Tutorial para Processamento de Dados de Single Cell RNA Sequencing da HIVE

Neste tutorial, vamos passar pelos passos necessários para processar dados de sequenciamento de RNA de célula única (scRNA-seq) provenientes do kit HIVE usando o softwear BeeNet. O objetivo é transformar arquivos FASTQ brutos em uma matriz de expressão gênica pronta para análise. Este tutorial é baseado no sistema Linux, como o Ubuntu, e portanto assume que você já tem esse sistema instalado. Caso não tenha, siga essas instruções.

0. Entendendo a estrutura dos dados

O primeiro passo para processar qualquer tipo de dados é entender sua estrutura. Quando fazemos um sequenciamento de scRNA, cada amostra pode vir dividida até 16 arquivos por amostra, sendo 2 arquivos de reads para até 8 lanes.

Durante o processo de sequenciamento, uma amostra pode ser distribuída em até 8 fileiras do equipamento, chamadas de lanes. Cada lane contém uma parte da amostra, dividida para facilitar o processamento simultâneo. Por isso, arquivos com a nomenclatura L001 - L008 representam a mesma amostra, dividida em até 8 lanes.

Além disso, cada lane gera dois arquivos de reads: R1 e R2, que contêm informações complementares sobre a amostra. Cada lane, portanto, possui um par de arquivos R1/R2

Assim, seus arquivos terão nomes no formato: Sample_L001_R1.fastq.gz, Sample_L001_R2.fastq.gz, Sample_L002_R1.fastq.gz, Sample_L001_R2.fastq.gz e assim por diante.

Independentemente do lane ou read, todos os arquivos FASTQ seguem a mesma estrutura. Cada read é composto por quatro linhas:

@Identificador único (nome da leitura)
Sequência (sequência de nucleotídeos)
+ (separador)
Score de qualidade (qualidade de cada base da sequência)

Nos arquivos R1 a sequência corresponde ao código de barras da célula e o UMI (Identificadores Moleculares Únicos) do mRNA, que são usados para identificar a célula de origem e quantificar as moléculas de mRNA. Já nos arquivos R2, a sequência corresponde ao mRNA transcrito da célula.

Exemplo de um read do arquivo R1:

@LH00136:22:22NGJMLT3:1:1106:48285:1000 1:N:0:AATCCGTT+AGCCTATT
CNCCGGAGCTGTTCCCTCGTTCTTTT
+
I#IIIIIIIIIIIIIIIIIIIIIIII

Exemplo de um read do arquivo R2:

@LH00136:22:22NGJMLT3:1:1106:48285:1000 2:N:0:AATCCGTT+AGCCTATT
ANCTCCGTGGATNCCTTCTNTTTGAATTTTTTTCTATAGNCTTGGTACTGT
+
I#I9-9IIIIII#IIIIII#IIIIIIIIIIIIIIIIIII#IIIIIIIIIII

Para ver as N primeiras linhas do seu arquivo fastq.gz, utilize o seguinte comando:

zcat caminho/para/arquivo.fastq.gz | head -n N

1. Configurando o Ambiente

Baixe os arquivos necessários para a construção do genoma referência: Acesse o site de Current Release do PlasmoDB, baixe os arquivos fasta e gff do genoma a ser utilizado.
Baixe e instale o BeeNet: Após se registrar no site da Honeycomb você receberá um link no seu email para baixar o software BeeNet. Copie o link e execute

wget (link)
chmod +x beenet

Isso dará permissão para o BeeNet atuar na pasta em que se encontra. Lembre-se que o BeeNet não funcionará em outra pasta

2. Convertendo GFF para GTF (se necessário)

Instale o gffread:

sudo apt-get install gffread

Converta GFF para GTF:

gffread your_annotation.gff -T -o your_annotation.gtf

Modifique GTF para BeeNet: o BeeNet exige que os arquivos gtf tenham \(;\) no final de todas as linhas, para isso rode o seguinte código em Python

# Python script para adicionar ';' ao final das linhas se não tiver

def add_semicolon_to_gtf(input_file, output_file):
    with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
        for line in infile:
            # Ignora linhas de comentário
            if line.startswith('#'):
                outfile.write(line)
            else:
                # Checa se a linha termina com (;) 
                if not line.strip().endswith(';'):
                    # Adiciona (;) se necessário
                    outfile.write(line.strip() + ';\n')
                else:
                    # Mantém a linha como está se já tiver (;)
                    outfile.write(line)

# Usar a função
input_gtf = 'your_input.gtf'   # Substitua por seu arquivo gtf
output_gtf = 'your_output.gtf' # Substitua pelo nome desejado para o novo arquivo gtf

add_semicolon_to_gtf(input_gtf, output_gtf)

print("Processing complete. New GTF file saved as", output_gtf)

Você pode também executar essa linha de comando

awk '{if ($0 !~ /;$/) $0=$0";"}1' seu_arquivo.gtf > seu_arquivo_modificado.gtf

3. Construindo o Genoma de Referência

Crie o genoma de referência com o beenet make-ref:

./beenet make-ref --ref=/caminho/para/genoma.fasta \
                 --gtf=/caminho/para/genoma.gtf \
                 --out=/caminho/para/pasta/output

4. Processando Arquivos FASTQ (opcional)

É muito provável que seus dados venham dividos em múltiplas lanes. O BeeNet lida com isso de forma automática, mas se você quiser juntar todas as lanes em um único arquivo, use o comando cat:

cat L1_R1.fastq.gz L2_R1.fastq.gz ... L8_R1.fastq.gz > sample_R1.fastq.gz
cat L1_R2.fastq.gz L2_R2.fastq.gz ... L8_R2.fastq.gz > sample_R2.fastq.gz

Você pode ver as primeiras linhas de cada arquivo utilizando

zcat sample_R1.fastq.gz | head

5. Criando arquivo txt

Ao invés de você passar os caminhos para os arquivos fastq na linha comando, o BeeNet permite que você crie um arquivo txt com esses caminhos. O arquivo deve conter um par de lane por linha (R1 + R2) e pode conter uma lane por linha (por isso não é necessário juntar as lanes como descrito no passo 4). Um exemplo desse arquivo é

Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L001_R1_001.fastq.gz Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L001_R2_001.fastq.gz
Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L002_R1_001.fastq.gz Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L002_R2_001.fastq.gz

Note que os arquivos da mesma lane se encontram na mesma linha separados por um espaço. Os caminhos devem ser relativos à pasta onde se está trabalhando

6. Realizando o Controle de Qualidade (QC)

Antes de executar o alinhamento e análise é essencial checar a qualidade dos seus arquivos FASTQ para se assegurar que você está trabalhando com bons dados.

Instale FastQC:

sudo apt-get install fastqc

Execute FastQC:

fastqc /path/to/R1.fastq.gz
fastqc /path/to/R2.fastq.gz

Revise o relatório de qualidade: Uma vez que o FastQC terminar, ele gerará um relatório HTML para cada arquivo. Abra esses relatórios no navegador e foque em:
- Qualidade por base da sequência: Procure por qualidade geral alta, especialmente no R1 (barcode e UMI).
- Conteúdo de adaptadores: Verifique se há contaminação.
- Conteúdo GC por sequência: Deve ser aproximadamente uniforme em R1 e pode conter um leve desvio em R2, visto que Plasmodium possui um leve viés para sequências AT.
- Nível de Sequências Duplicadas: É normal em scRNA-seq ter sequências duplicadas.

7. Executando o BeeNet

Execute beenet analyze para alinhar as leituras e gerar a matriz gene-célula:

./beenet analyze --sample-name=nome/para/output \
                 --num-barcodes=num \
                 --ref=caminho/para/genoma/referencia \
                 --fastqs=caminho/para/arquivo/txt \
                 --out=nome/pasta/output

--sample-name é o nome do arquivo que será criado no output

--num-barcodes é o número esperado de barcodes na amostra (aproximadamente 40% das células de input)

--ref é o caminho para a pasta criada pelo output de make-ref

--fastqs é o caminho para o arquivo txt criado com os nomes dos arquivos fastq

--out é o nome da pasta de output

Importante: ao executar analyze você pode ou usar o argumento –fastqs com o arquivo txt ou passar os caminhos no argumento –ref. Exemplo:

./beenet analyze --sample-name=Beto_P_knowlesi_cultura_80K_1 --num-barcodes=32000 --ref=Genomes/Pk-68/PlasmoDB-68_PknowlesiH --fastqs=Pk1/fastq/fastq.txt --out=Pknowlesi_1_folder

# OU #

./beenet analyze --sample-name=Beto_P_knowlesi_cultura_80K_1 --num-barcodes=32000 --ref=Genomes/Pk-68/PlasmoDB-68_PknowlesiH Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L001_R1_001.fastq.gz Pk1/fastq/Beto_P_knowlesi_cultura_80K_1_S1_L001_R2_001.fastq.gz --out=Pknolesi_1_folder

8. Entendendo o output

Seu output será uma pasta contendo os seguintes arquivos:

bam: usado para criar as count matrix
RCM: Read Count Matrix
TCM: Transcript Count Matrix
CMSummary: Count Matrix Summary
ReadsQC: Controle de Qualidade das reads
SampleQC: Controle de Qualidade da amostra
log: log da execução

Você deve estar se perguntando qual é a diferença entre RCM e TCM. O arquivo RCM contém a matriz de contagem bruta de leituras mapeadas para cada gene em cada célula. Esse arquivo reflete o número total de leituras antes da normalização ou deduplicação, incluindo leituras duplicadas causadas pela amplificação do material genético. Já o arquivo TCM contém a matriz de contagem de transcritos, onde as leituras são deduplicadas com base nos UMIs, fornecendo uma estimativa mais precisa da abundância de transcritos únicos por gene em cada célula. Ou seja, enquanto o RCM representa os dados brutos diretamente do sequenciamento, o TCM é utilizado nas análises biológicas por refletir de forma mais precisa a expressão gênica real após a remoção de duplicatas técnicas.

9. Pós-processamento no R

Após a execução do beenet analyze, você pode carregar os resultados no R para análise posterior usando pacotes como Seurat.

O script a seguir funciona se você tiver somente uma amostra, uma amostra dividida em mais que um dataset e também se você tiver diferentes amostras biológicas que quer analisar em conjunto (como um mesmo objeto). Esse último caso pode ser verdade para dados de diferentes estágios do parasita que você pretende analisar como se fosse um conjunto de dados único.

1. Configurando o ambiente e metadados

O primeiro passo é criar uma lista com os nomes das amostras (aqui amostra significa o set de dados gerados). Os nomes das amostras nesta lista devem corresponder aos nomes das amostras nos arquivos de saída do BeeNet. Por exemplo, se o arquivo da Matriz de Contagem de Transcrição (TCM) for nomeado como Amostra1_20230426182713_TCM.tsv.gz, então o nome da sua amostra é Amostra1.

Na lista groupslist você pode adicionar qualquer metadado que queira incorporar à amostra. Se, por exemplo, a Amostra1 corresponde a uma amostra biológica do tipo selvagem e a Amostra2 corresponde a uma amostra biológica com tratamento, você pode colocar ('wild_type', 'treatment') no groupslist. Independente do que você colocar nessa lista, será tratado como mais um metadado da amostra em questão.

Só lembre-se de se certificar que a ordem das amostras e dos metadados está correta!

# Set working directory to folder where your files are
setwd('caminho/para/pasta/com/arquivos')

# Transcript count matrix files
CMfile <- list.files(path = c("./Amostra1", "./Amostra2", "./Amostra3"), pattern = "TCM.tsv.gz", full.names = TRUE)

# ReadsQC files
readsFile <- list.files(path = c("./Amostra1", "./Amostra2", "./Amostra3"), pattern = "ReadsQC.tsv", full.names = TRUE)

# SampleQC file
samplesFile <- list.files(path = c("./Amostra1", "./Amostra2", "./Amostra3"), pattern = "SampleQC.tsv", full.names = TRUE)

# Certifique-se que a ordem das amostras e dos grupos está correta!

sampleslist <- c("Amostra1", "Amostra2", "Amostra3")
groupslist <- c("wild_type", "treatment_1", "treatment_2")

Agora você deve ver uma lista das suas amostras e seus grupos no ambiente do R.

2. Lendo os dados e criando um objeto Seurat

Aqui definimos uma função chamada create_S.O() para criar objetos Seurat a partir das saídas BeeNet. Esta função lê matrizes de contagem e adiciona metadados ao seu objeto Seurat.

create_S.O<- function(sample, group_name = ""){

  # Carregar matriz TCM
  dataTran <- read.table(CMfile[grep(sample, CMfile)],
                         header=1,
                         row.names=1,
                         check.names=FALSE , 
                         sep="\t", 
                         quote = "", 
                         stringsAsFactors=0) 
  
  # Carregar informação dos reads de cada código de barras
  dataReads <- read.table(readsFile[grep(sample, readsFile)],
                        header=1,
                        row.names=1)
  
  # Certifique-se que a ordem das informações de reads é a mesma que a ordem da matriz TCM
  table(rownames(dataReads) == colnames(dataTran))
  reads=dataReads[rownames(dataReads) %in% colnames(dataTran),]
  
  # Criar a Seurat Object
  # Aqui min.cells significa o número mínimo de células que devem ter um determinado gene para que ele seja válido e min.features é o número mínimo de genes em uma célula para que ela seja válida
  obj <- CreateSeuratObject(counts=dataTran, project = sample, min.cells = 10, min.features = 100)
  
  # Calcular metadados
  logtotreads <- log10(as.matrix(reads)[,1])
  readAll <- dataReads[,1]
  readMap <- reads[,2]
  readMapPct <- readMap/readAll
  readUnmap <- readAll-readMap
  readUnmapPct <- readUnmap/readAll
  readExon <- reads[,3]
  readExonPct <- readExon/readAll
  readExonMapPct <- readExon/readMap
  names(readExon) <- colnames(obj)
  obj <- AddMetaData(obj, readExon, "ExonReads")
  names(readExonMapPct) <- colnames(obj)
  names(readExonPct) <- colnames(obj)
  names(readMap) <- colnames(obj)
  names(readAll) <- colnames(obj)
  
  # Adicionar metadados ao objeto
  obj <- AddMetaData(obj, logtotreads, "Log.TotReads") 
  obj <- AddMetaData(obj, readExonMapPct, "ExonvMapped")
  obj <- AddMetaData(obj, readExonPct, "ExonvTotal")
  obj <- AddMetaData(obj, readMap, "reads.mapped")
  obj <- AddMetaData(obj, readAll, "reads.Total")

  # Calcular complexidade e adicionar ao objeto
  comp=obj@meta.data$ExonReads/obj@meta.data$nCount_RNA
  names(comp)<-rownames(obj@meta.data)
  obj <- AddMetaData(obj, comp, "Complexity")
  
  # Adicionar grupo como metadado (aqueles que você colocou em groupslist)
  obj$group <- group_name
  obj$groupSamp <- paste(group_name, sample, sep = "_")

  return(obj)
}

Em seguida, execute essa função em cada uma das suas amostras para gerar um objeto Seurat para cada. A partir daí, se necessário, combine os objetos individuais em objetos finais.

# Crie objetos Seurat para cada amostra
# Essa etapa pode demorar um pouco se você tiver muitas amostras
seurat_list <- purrr::map2(sampleslist, groupslist, 
                    ~ create_S.O(sample = .x,
                                             group_name = .y))

# Se você tiver mais que uma amostra que deseja combinar em um único objeto Seurat
if(length(seurat_list) > 1) {
  S.O.RNA <- merge(x = seurat_list[[1]], 
             y = seurat_list[2:length(seurat_list)], 
             add.cell.ids = sampleslist)
} else {
  S.O.RNA <- seurat_list[[1]]
}

# Caso você tenha combinado as amostras, remova a lista de objetos individuais para liberar a memória
remove(seurat_list)

3. Limpando e pré-processando os dados

Agora que você possui seu objeto Seurat, chamado S.O.RNA, o próximo passo é realizar a limpeza e normalização dos dados. A limpeza envolve a remoção de células com valores extremos em relação à quantidade de genes expressos (nFeature) e ao número total de UMIs (nCount). Isso é importante, pois células com valores muito altos provavelmente são duplicatas, enquanto células com valores muito baixos podem ser indicativas de células mortas ou de baixa qualidade. Para esse processo, utilizamos limites baseados nos percentis de 10% e 98%. Esses percentis são escolhidos porque representam valores que excluem as células mais extremas sem afetar a maioria dos dados.

VlnPlot(S.O.RNA, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
FeatureScatter(S.O.RNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
quant1 <- quantile(S.O.RNA$nFeature_RNA, probs = c(0.1, 0.98))
quant2 <- quantile(S.O.RNA$nCount_RNA, probs = c(0.1, 0.98))
S.O.RNA <- subset(S.O.RNA, subset = nFeature_RNA > quant1[1] & nFeature_RNA < quant1[2] & 
                    nCount_RNA > quant2[1] & nCount_RNA < quant2[2])
table(S.O.RNA$orig.ident)

A função prep_S.O() realiza diversas etapas importantes no processamento dos dados: normaliza os dados, identifica genes outliers, escala e centraliza os valores, executa uma análise de componentes principais (PCA), encontra células semelhantes (vizinhas), as clusteriza e, por fim, realiza uma análise UMAP para visualização.

prep_S.O <- function(S.O, res = 0.1, var.features = F, down.sample = F){
  set.seed(100)
  S.O <- NormalizeData(S.O, normalization.method = "LogNormalize", scale.factor = 10000)
  S.O <- FindVariableFeatures(S.O, selection.method = "vst", nfeatures = 6000)
  if(var.features){
    ## Work on variable features only
    S.O <- subset(S.O, features = VariableFeatures(S.O))
  }
  all.genes <- rownames(S.O)
  S.O <- ScaleData(S.O, features = all.genes)
  S.O <- RunPCA(S.O, features = VariableFeatures(object = S.O))
  S.O <- FindNeighbors(S.O, dims = 1:10, reduction = 'pca')
  S.O <- FindClusters(S.O, resolution = res)
  if(down.sample){
    S.O <- subset(x = S.O, downsample = 800)
    S.O <- FindNeighbors(S.O, dims = 1:13)
    S.O <- FindClusters(S.O, resolution = res)
  }
  S.O <- RunUMAP(S.O, dims = 1:13, n.components = 3L)
  return(S.O)
}

S.O.RNA <- prep_S.O(S.O.RNA)

# Plotar gráficos de PCA e UMAP
DimPlot(S.O.RNA, reduction = 'pca')
DimPlot(object = S.O.RNA, label = TRUE, reduction = 'umap') + NoLegend()

# Salve seu objeto para referência futura!
saveRDS(S.O.RNA, '../S.O.exp.rds')

4. Passo final

Agora que seu objeto está limpo, normalizado e pré-processado, ele está pronto para você fazer todas as análises que quiser! Divirta-se!