DGE and DTE analysis of Salmon/Kallisto inputs using Sleuth
Analysis and result presented was performed with Salmon counts, Code snippet to import Kallisto counts is also provided
DGE analysis
cd/home/USER/SSAPs/analysis/dgeR
Load required libraries
library(data.table)library(sleuth)
Import clinical data
Load clinical information and define file path
dir ="../../"sampleData =paste0(dir, "clinical.txt")sampleData =fread(sampleData)rownames(sampleData) = sampleData$ENA_RUNsampleData$individual =as.factor(sampleData$individual)sampleData$paris_classification =as.factor(sampleData$paris_classification)# Use relevel() to set adjacent normal samples as referencesampleData$paris_classification =relevel(sampleData$paris_classification, ref ="normal")
Construct t2g data.frame
Load the annoData data.frame that contain per-transcript annotation derived from GENCODE GTF file and constrct a transcript-to-gene data.frame t2g in the specified format by Sleuth
> s2c sample condition pathERR2675454 ERR2675454 0-IIa ../../salmon-bs/ERR2675454ERR2675455 ERR2675455 normal ../../salmon-bs/ERR2675455ERR2675458 ERR2675458 0-IIa ../../salmon-bs/ERR2675458ERR2675459 ERR2675459 normal ../../salmon-bs/ERR2675459ERR2675460 ERR2675460 0-IIa ../../salmon-bs/ERR2675460ERR2675461 ERR2675461 normal ../../salmon-bs/ERR2675461ERR2675464 ERR2675464 0-IIa ../../salmon-bs/ERR2675464ERR2675465 ERR2675465 normal ../../salmon-bs/ERR2675465ERR2675468 ERR2675468 0-IIa ../../salmon-bs/ERR2675468ERR2675469 ERR2675469 normal ../../salmon-bs/ERR2675469ERR2675472 ERR2675472 0-IIa ../../salmon-bs/ERR2675472ERR2675473 ERR2675473 normal ../../salmon-bs/ERR2675473ERR2675476 ERR2675476 0-IIa ../../salmon-bs/ERR2675476ERR2675477 ERR2675477 normal ../../salmon-bs/ERR2675477ERR2675478 ERR2675478 0-IIa ../../salmon-bs/ERR2675478ERR2675479 ERR2675479 normal ../../salmon-bs/ERR2675479ERR2675480 ERR2675480 0-IIa ../../salmon-bs/ERR2675480ERR2675481 ERR2675481 normal ../../salmon-bs/ERR2675481ERR2675484 ERR2675484 0-IIa ../../salmon-bs/ERR2675484ERR2675485 ERR2675485 normal ../../salmon-bs/ERR2675485
Sleuth will compare all conditions against the reference level (i.e. first level in the factor variable). Use the levels() function to check the order of the condition and userelevel() function to reassign reference level if necessary
> s2c sample condition pathERR2675454 ERR2675454 0-IIa ../../kallisto/ERR2675454ERR2675455 ERR2675455 normal ../../kallisto/ERR2675455ERR2675458 ERR2675458 0-IIa ../../kallisto/ERR2675458ERR2675459 ERR2675459 normal ../../kallisto/ERR2675459ERR2675460 ERR2675460 0-IIa ../../kallisto/ERR2675460ERR2675461 ERR2675461 normal ../../kallisto/ERR2675461ERR2675464 ERR2675464 0-IIa ../../kallisto/ERR2675464ERR2675465 ERR2675465 normal ../../kallisto/ERR2675465ERR2675468 ERR2675468 0-IIa ../../kallisto/ERR2675468ERR2675469 ERR2675469 normal ../../kallisto/ERR2675469ERR2675472 ERR2675472 0-IIa ../../kallisto/ERR2675472ERR2675473 ERR2675473 normal ../../kallisto/ERR2675473ERR2675476 ERR2675476 0-IIa ../../kallisto/ERR2675476ERR2675477 ERR2675477 normal ../../kallisto/ERR2675477ERR2675478 ERR2675478 0-IIa ../../kallisto/ERR2675478ERR2675479 ERR2675479 normal ../../kallisto/ERR2675479ERR2675480 ERR2675480 0-IIa ../../kallisto/ERR2675480ERR2675481 ERR2675481 normal ../../kallisto/ERR2675481ERR2675484 ERR2675484 0-IIa ../../kallisto/ERR2675484ERR2675485 ERR2675485 normal ../../kallisto/ERR2675485
Construct Sleuth object
Construct the Sleuth object so with sleuth_prep. This object will store the metadata, model to be used for differential testing, transcript-to-gene and bootstrap information
Use aggregation_column = "ens_gene" to allow summarize the data on the gene level and gene_mode = TRUE to do counts aggregation at the gene level for normalization, transformation, and downstream steps.
The default transformation of counts is natural log, i.e. log(x+0.5). We change this to log base 2 by assigning an appropriate transformation_function
reading in kallisto resultsdropping unused factor levels....................normalizing est_counts51524 targets passed the filternormalizing tpmmerging in metadataaggregating by column: ens_gene17692 genes passed the filtersummarizing bootstraps....................
Differential expression analysis
The likelihood ratio test is performed with a 3-step procedure. First, fit the full model by fit_name = "full", then performs a second fit to a “reduced” model that presumes abundances are equal in the two conditions. Finally, Sleuth will perform LRT to identify transcripts with a significantly better fit with the “full” model
We construct a normData data.frame to store per-group normalised mean transcripts per million (TPM) and normalised TPM of all samples, and a deData data.frame to store the Sleuth LRT results. We merge both data.frame with a annoData data.frame that contain per-transcript annotation derived from GENCODE GTF file before exporting the results in tabular format
> normData[1:10,1:14] GeneID GeneSymbol Chromosome Start End Class Strand Length normal 0-IIa ERR2675454 ERR2675455 ERR2675458 ERR2675459
24914 ENSG00000223972.5 DDX11L1 chr1 11869 14409 transcribed_unprocessed_pseudogene + 2540 0.02578694 0.002660157 0.0000000 0.0000000 0.000000 0.000000
27153 ENSG00000227232.5 WASH7P chr1 14404 29570 unprocessed_pseudogene - 15166 9.64654711 7.418128917 5.0318657 7.0457119 5.043944 11.669079
54425 ENSG00000278267.1 MIR6859-1 chr1 17369 17436 miRNA - 67 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
36601 ENSG00000243485.5 MIR1302-2HG chr1 29554 31109 lncRNA + 1555 0.01174571 0.000000000 0.0000000 0.1174571 0.000000 0.000000
57122 ENSG00000284332.1 MIR1302-2 chr1 30366 30503 miRNA + 137 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
34149 ENSG00000237613.2 FAM138A chr1 34554 36081 lncRNA - 1527 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
49448 ENSG00000268020.3 OR4G4P chr1 52473 53312 unprocessed_pseudogene + 839 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
35334 ENSG00000240361.2 OR4G11P chr1 57598 64116 transcribed_unprocessed_pseudogene + 6518 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
16136 ENSG00000186092.6 OR4F5 chr1 65419 71585 protein_coding + 6166 0.00000000 0.000000000 0.0000000 0.0000000 0.000000 0.000000
34420 ENSG00000238009.6 AL627309.1 chr1 89295 133723 lncRNA - 44428 0.37463259 0.356671874 0.1921258 0.5975709 0.389049 0.307692
> deData[1:10,] GeneID GeneSymbol Chromosome Start End Class Strand Length baseMean log2fc pvalue padj
24914 ENSG00000223972.5 DDX11L1 chr1 11869 14409 transcribed_unprocessed_pseudogene + 2540 NA -3.27705754 NA NA
27153 ENSG00000227232.5 WASH7P chr1 14404 29570 unprocessed_pseudogene - 15166 5.289557 -0.37895730 0.1989157 0.2871868
54425 ENSG00000278267.1 MIR6859-1 chr1 17369 17436 miRNA - 67 NA NA NA NA
36601 ENSG00000243485.5 MIR1302-2HG chr1 29554 31109 lncRNA + 1555 NA NA NA NA
57122 ENSG00000284332.1 MIR1302-2 chr1 30366 30503 miRNA + 137 NA NA NA NA
34149 ENSG00000237613.2 FAM138A chr1 34554 36081 lncRNA - 1527 NA NA NA NA
49448 ENSG00000268020.3 OR4G4P chr1 52473 53312 unprocessed_pseudogene + 839 NA NA NA NA
35334 ENSG00000240361.2 OR4G11P chr1 57598 64116 transcribed_unprocessed_pseudogene + 6518 NA NA NA NA
16136 ENSG00000186092.6 OR4F5 chr1 65419 71585 protein_coding + 6166 NA NA NA NA
34420 ENSG00000238009.6 AL627309.1 chr1 89295 133723 lncRNA - 44428 1.263687 -0.07087897 0.2654303 0.3487556
Exploring results
Principal component plot of the samples
png("DGE_sleuth_PCA.png", width=7, height=7, units ="in", res =300)plot_pca(so, color_by ="condition", text_labels =TRUE, units ="tpm")dev.off()
Plot sample heatmap using the Jensen-Shannon divergence
png("DGE_sleuth_sample_heatmap.png", width=7, height=7, units ="in", res =300)plot_sample_heatmap(so, use_filtered =FALSE)dev.off()
Lower divergence values represent samples that are more similar to each other
Plot clustered transcript heatmap
Create expression heatmap for top 40 genes with lowest q-values
png("DGE_sleuth_transcript_heatmap.png", width=7, height=7, units ="in", res =300)plot_transcript_heatmap(so, cluster_transcripts =TRUE,transcripts =subset(res[order(res$qval),], qval <0.05)$target_id[1:40]) # top40dev.off()
Plot the normalized bootstraps across all samples
Plot expression with bootstrap variation for selected gene
png("DGE_sleuth_plot_bootstrap.png", width=7, height=7, units ="in", res =300)plot_bootstrap(so, target_id ="ENSG00000062038.14", units ="tpm", color_by ="condition")dev.off()
dir ="../../"sampleData =paste0(dir, "clinical.txt")sampleData =fread(sampleData)rownames(sampleData) = sampleData$ENA_RUNsampleData$individual =as.factor(sampleData$individual)sampleData$paris_classification =as.factor(sampleData$paris_classification)# Use relevel() to set adjacent normal samples as referencesampleData$paris_classification =relevel(sampleData$paris_classification, ref ="normal")
Construct t2g data.frame
Load the annoData data.frame that contain per-transcript annotation derived from GENCODE GTF file and constrct a transcript-to-gene data.frame t2g in the specified format by Sleuth
> s2c sample condition pathERR2675454 ERR2675454 0-IIa ../../salmon-bs/ERR2675454ERR2675455 ERR2675455 normal ../../salmon-bs/ERR2675455ERR2675458 ERR2675458 0-IIa ../../salmon-bs/ERR2675458ERR2675459 ERR2675459 normal ../../salmon-bs/ERR2675459ERR2675460 ERR2675460 0-IIa ../../salmon-bs/ERR2675460ERR2675461 ERR2675461 normal ../../salmon-bs/ERR2675461ERR2675464 ERR2675464 0-IIa ../../salmon-bs/ERR2675464ERR2675465 ERR2675465 normal ../../salmon-bs/ERR2675465ERR2675468 ERR2675468 0-IIa ../../salmon-bs/ERR2675468ERR2675469 ERR2675469 normal ../../salmon-bs/ERR2675469ERR2675472 ERR2675472 0-IIa ../../salmon-bs/ERR2675472ERR2675473 ERR2675473 normal ../../salmon-bs/ERR2675473ERR2675476 ERR2675476 0-IIa ../../salmon-bs/ERR2675476ERR2675477 ERR2675477 normal ../../salmon-bs/ERR2675477ERR2675478 ERR2675478 0-IIa ../../salmon-bs/ERR2675478ERR2675479 ERR2675479 normal ../../salmon-bs/ERR2675479ERR2675480 ERR2675480 0-IIa ../../salmon-bs/ERR2675480ERR2675481 ERR2675481 normal ../../salmon-bs/ERR2675481ERR2675484 ERR2675484 0-IIa ../../salmon-bs/ERR2675484ERR2675485 ERR2675485 normal ../../salmon-bs/ERR2675485
Sleuth will compare all conditions against the reference level (i.e. first level in the factor variable). Use the levels() function to check the order of the condition and userelevel() function to reassign reference level if necessary
> s2c sample condition pathERR2675454 ERR2675454 0-IIa ../../kallisto/ERR2675454ERR2675455 ERR2675455 normal ../../kallisto/ERR2675455ERR2675458 ERR2675458 0-IIa ../../kallisto/ERR2675458ERR2675459 ERR2675459 normal ../../kallisto/ERR2675459ERR2675460 ERR2675460 0-IIa ../../kallisto/ERR2675460ERR2675461 ERR2675461 normal ../../kallisto/ERR2675461ERR2675464 ERR2675464 0-IIa ../../kallisto/ERR2675464ERR2675465 ERR2675465 normal ../../kallisto/ERR2675465ERR2675468 ERR2675468 0-IIa ../../kallisto/ERR2675468ERR2675469 ERR2675469 normal ../../kallisto/ERR2675469ERR2675472 ERR2675472 0-IIa ../../kallisto/ERR2675472ERR2675473 ERR2675473 normal ../../kallisto/ERR2675473ERR2675476 ERR2675476 0-IIa ../../kallisto/ERR2675476ERR2675477 ERR2675477 normal ../../kallisto/ERR2675477ERR2675478 ERR2675478 0-IIa ../../kallisto/ERR2675478ERR2675479 ERR2675479 normal ../../kallisto/ERR2675479ERR2675480 ERR2675480 0-IIa ../../kallisto/ERR2675480ERR2675481 ERR2675481 normal ../../kallisto/ERR2675481ERR2675484 ERR2675484 0-IIa ../../kallisto/ERR2675484ERR2675485 ERR2675485 normal ../../kallisto/ERR2675485
Construct Sleuth object
Construct the Sleuth object so with sleuth_prep. This object will store the metadata, model to be used for differential testing, transcript-to-gene and bootstrap information.
The default transformation of counts is natural log, i.e. log(x+0.5). We change this to log base 2 by assigning an appropriate transformation_function
reading in kallisto resultsdropping unused factor levels....................normalizing est_counts51524 targets passed the filternormalizing tpmmerging in metadatasummarizing bootstraps....................
Differential expression analysis
The likelihood ratio test is performed with a 3-step procedure. First, fit the full model by fit_name = "full", then performs a second fit to a “reduced” model that presumes abundances are equal in the two conditions. Finally, Sleuth will perform LRT to identify transcripts with a significantly better fit with the “full” model
We construct a normData data.frame to store per-group normalised mean transcripts per million (TPM) and normalised TPM of all samples, and a deData data.frame to store the Sleuth LRT results. We merge both data.frame with a annoData data.frame that contain per-transcript annotation derived from GENCODE GTF file before exporting the results in tabular format