Worklfow anàlisis RNA-seq
- Aquest document descriu el worflow per analitzar les dades de RNA-seq, des de el seu alineament fins als anàlisis d’enriquiment funcional.
- Conté codi BASH, identificat en el codi amb el
chunk {bash} i en l’informe amb color verd pàlid, i codi
R, identificar en el codi amb el chunk {r} i
amb color violeta pàlid en l’informe. Els programes amb BASH no
s’executen a R, van directament a la Terminal/Consola. Son els segünets:
- FastQC
- MultiQC
- STAR
- samtools
- El codi esta estructurat dins chunks, si es vol executar una part del codi cal substituïr l’atribut eval=F per eval=T.
- Està preparat per ser executat en l’ordenador de Genòmica a través del link facilitat. En cas de que es vulgui instalar en un altra ordenador caldrà instalar les dependències corresponents.
Quality Control of FASTQ files
El primer pas del flux de treball RNA-Seq és agafar els fitxers FASTQ rebuts de la seqüenciació i avaluar la qualitat de les lectures de la seqüència.
Assessing quality with FastQC
mkdir ./RESULTATS/QC
fastqc -o ./RESULTATS/QC ./RAW/*/*.gz# Es borrarà els QC anteriors
rm -d ./RESULTATS/QC/multiqc_data/*
rm -d ./RESULTATS/QC/multiqc_data/
rm ./RESULTATS/QC/multiqc_report.html
multiqc -o ./RESULTATS/QC/ ./RESULTATS/QC/ Interpretció FastQC Report
Basic statistics
Estadístiques bàsiques de la mostra. Comprovar que la longitud de lectura i el contingut %GC siguin els esperats.
Per base sequence quality
Gràfica amb la distribució de les puntuacions de qualitat a cada posició de la lectura en totes les lectures. L’eix Y dóna les puntuacions de qualitat, mentre que l’eix X representa la posició a la lectura. La codificació de colors de la trama indica la qualitat de la puntuació (Alt-Mitg-Baix).
Pot indicar possibles errors técnicns durant la seqüenciació.
Per sequence quality scores
Mostra la puntuació de qualitat mitjana a l’eix X i el nombre de seqüències amb aquesta mitjana a l’eix Y. S’espera que la majoria de les lectures tinguin una puntuació de qualitat mitjana alta.
Per base sequence content
Sempre dóna un FAIL per a les dades de la seqüència d’ARN. Això es deu al fet que les primeres 10-12 bases resulten de l’encebament de l’hexàmer “atzar” que es produeix durant la preparació de la biblioteca d’ARN-seq. Aquest cebament no és tan aleatori com podríem esperar, donant un enriquiment en bases particulars per a aquests nucleòtids inicials.
Per sequence GC content
Mostra la distribució de GC en totes les seqüències. La distribució hauria de ser normal tret que hi hagi seqüències sobrerepresentades (pics aguts en una distribució normal) o contaminació amb un altre organisme (pic ampli). Les seqüències sobrerepresentades poden indicar contaminació o un gen molt sobreexpressat.
Over-represented sequences
Presència de seqüències (almenys 20 pb) que es produeixen en més del 0,1% del nombre total de seqüències. Aquesta taula ajuda a identificar la contaminació, com ara seqüències de vectors o adaptadors. Si el contingut %GC falla, aquesta taula pot ajudar a identificar la font. Si no apareix com a adaptador o vector conegut, pot ajudar a BLAST la seqüència per determinar la identitat.
Adapter Content
Avís si hi ha alguna seqüència en més del 10% de totes les lectures.
FastQC Report de totes les mostres (MultiQC)
`
Alineament
Index
STAR Aligner
Creating a genome index
20 minuts
Genoma de referència i GTF descarregats de http://www.ensembl.org/info/data/ftp/index.html
- Mus_musculus.GRCm39.dna.primary_assembly.fa
- Mus_musculus.GRCm39.105.gtf
–readFilesCommand zcat: No funciona en el genoma de referenica i el gtf
–sjdbOverhang option needs to be specified for detecting possible splicing sites. Default default: 100 Hauria de ser longitud de lectura -1
S’ha de fer correr una única vegada. Un cop creat el index (50 minuts aprox) evitar tornar-lo a executar (eval=F).
–runThreadN 10: No tocar
STAR --runMode genomeGenerate --genomeDir REF/ --genomeFastaFiles REF/Mus_musculus.GRCm39.dna.primary_assembly.fa --sjdbGTFfile REF/Mus_musculus.GRCm39.105.gtf --outFileNamePrefix RESULTATS/INDEX/INDEX --runThreadN 10 --sjdbOverhang 99Aligning reads to the genome
4~6 minuts per mostra.
Alineament mostra per mostra:
STAR --genomeDir RESULTATS/INDEX/ --readFilesIn RAW/C4E2/C4E2_1.fq.gz RAW/C4E2/C4E2_2.fq.gz --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix alignments/C4E2 --runThreadN 10
STAR --genomeDir RESULTATS/INDEX/ --readFilesIn RAW/C4E4/C4E4_1.fq.gz RAW/C4E4/C4E4_2.fq.gz --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix alignments/C4E4 --runThreadN 10
STAR --genomeDir RESULTATS/INDEX/ --readFilesIn RAW/C5E1/C5E1_1.fq.gz RAW/C5E1/C5E1_2.fq.gz --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix alignments/C5E1 --runThreadN 10
STAR --genomeDir RESULTATS/INDEX/ --readFilesIn RAW/C5E2/C5E2_1.fq.gz RAW/C5E2/C5E2_2.fq.gz --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix alignments/C5E2 --runThreadN 10Loop per alinear tores les mostres dins la carpeta RAW.
BAM to SAM
No cal. STAR ofereix la sortida directament a BAM (si es vol). Si es necessita obtenir el SAM per algún motiu es pot fer amb la segünet instrucció (canviar eval=F per eval=T)
samtools view -h alignments/TEST_5Aligned.sortedByCoord.out.bam > alignments/TEST_5.samVisualització IGV, browser etc
Fer l’index de cada mostra. Es crea C4E2Aligned.sortedByCoord.out.bam.bai
samtools index alignments/C4E2Aligned.sortedByCoord.out.bam
Anàlisi qualitat dels alignment
10 min aprox per mostra
WARNING: out of memory!
Ull viu amb la memòria.
Aquest pas requereix bastanta memòria. Hi ha una carpeta tempora que després es buida sola. El resultat final es desa en la carpeta RESULTATS/QC
mkdir ../tmpexport _JAVA_OPTIONS="-Djava.io.tmpdir=../tmp -Xmx6G"
../qualimap_v2.2.1/qualimap rnaseq -pe -bam alignments/C4E4Aligned.sortedByCoord.out.bam -gtf REF/Mus_musculus.GRCm39.105.gtf -outdir RESULTATS/QC/C4E4Aligned -p strand-specific-reverse
# 5.1 GBMultiQC
# Aquest multiqc porta el nom multiqc_report_1.html.
# Es borrarà els QC anteriors
rm -d ./RESULTATS/QC/multiqc_data_1/*
rm -d ./RESULTATS/QC/multiqc_data_1/
rm ./RESULTATS/QC/multiqc_report_1.html
multiqc -o ./RESULTATS/QC/ ./RESULTATS/QC/ Expressió diferencial
# Lectura dels Counts
arxius_reads<-list.files("./alignments/",pattern = "ReadsPerGene.out.tab",full.names = T)
count_matrix = lapply(arxius_reads, function(x) {read.delim(x,header = F,col.names = c("gene ID","counts unstranded","counts 1st read strand","counts2nd read strand"),skip=4)})
columna_ID = colnames(head(count_matrix[[1]]))[c(4)]
result = lapply(count_matrix, "[",, columna_ID)
result<-(data.frame(result))
# Assignar els gens
rownames(result)<-count_matrix[[1]]$gene.ID
## reverse stranded -> 4ta columna# Crear taula fenotipica
sample_sheet<-
data.frame(
'filename'=basename(arxius_reads),
'sample name'=gsub("ReadsPerGene.out.tab","",basename(arxius_reads)),
'group'=c(1,1,2,2))
rownames(sample_sheet)<-sample_sheet$filename
# Guardar i modificar en l'excel
write_xlsx(sample_sheet,"./sample_sheet.xlsx")
# Lerctura de la taula modificada
sample_sheet<-read_excel("./sample_sheet.xlsx")
sample_sheet$group<-as.factor(sample_sheet$group)
datatable(sample_sheet,caption = "Taula fenotip per completar",rownames = F,extensions = 'Buttons',options = list(
pageLength = 5,dom = 'Bfrtip',buttons = c('copy', 'csv', 'excel')))colnames(result)<-sample_sheet$filename
# crear el DESeqDataSet. A Design s'ha de definir la comparativa
se_star_matrix <- DESeqDataSetFromMatrix(countData = result,
colData = sample_sheet,
design = ~ group)Filtratge
- Filtrar gens amb més de 10 reads
pre_n<-nrow(se_star_matrix)
# Filtrar gens amb més de 10 reads
se_star_matrix <- se_star_matrix[rowSums(counts(se_star_matrix)) > 10, ]
# Number of genes left after low-count filtering:
post_n<-nrow(se_star_matrix)
print(paste0("S'han eliminat ",pre_n-post_n, " gens"))[1] “S’han eliminat 37316 gens”
Model estadístic
Transformació de les dades
vsd <- vst(se_star_matrix_2, blind = FALSE)
rld <- rlog(se_star_matrix_2, blind = FALSE)
dds <- estimateSizeFactors(se_star_matrix_2)
df <- bind_rows(
as_data_frame(log2(counts(se_star_matrix_2, normalized=TRUE)[, 1:2]+1)) %>%
mutate(transformation = "log2(x + 1)"),
as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))
colnames(df)[1:2] <- c("x", "y")
lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
coord_fixed() + facet_grid( . ~ transformation) Visualització
resultsNames(se_star_matrix_2)## [1] "Intercept" "group_2_vs_1"
de <- results(object = se_star_matrix_2,
name=resultsNames(se_star_matrix_2)[2])
de_shrink <- lfcShrink(dds = se_star_matrix_2,coef=resultsNames(se_star_matrix_2)[2],type="apeglm")Gens DE significatius
dir.create("./RESULTATS/TAULES")
de_symbols<-merge(log_counts,unique(gtf_df[,c("gene_id","gene_name")]),by.x="row.names",by.y="gene_id")
de_symbols <- merge(unique(gtf_df[,c("gene_id","gene_name")]), data.frame(ID=rownames(de_shrink), de_shrink), by=1, all=F)
write.table(de_symbols, "./RESULTATS/TAULES/deseq2_results.txt", quote=F, col.names=T, row.names=F, sep="\t")
datatable(de_symbols %>% filter(padj<=0.05),caption = "Gene Ontology (ORA)",extensions = 'Buttons',options = list(
columnDefs = list(list(className = 'dt-center', targets = 5)),
pageLength = 5,dom = 'Bfrtip',buttons = c('copy', 'csv', 'excel')))Significància biologica. EXEMPLES NO SIGNIFICATIUS!!!!!
Rich Factor: proporció de gens introduïts (per exemple, DE) que s’assignen a un terme en relació a tots els gens estudiats anotats en aquest terme.
Gene Ratio: proporció de gens introduïts (per exemple, DE) que s’assignen a un terme en relació a tots els gens en aquest terme.
#Referències: