Hello everyone!

In my PhD thesis, I am working with differential expression analysis, trying to verify genes related to thermoregulation in beef cattle.
I have 66 sequenced, paired-end RNA samples, divided into 2 lanes per sample.
The first steps were to check the quality of these samples by FASTQC, then do the trimming (using the Trimmomatic software with the parameters: LEADING: 30 TRAILING: 30 SLIDINGWINDOW: 5: 20 AVGQUAL: 20 MINLEN: 50) and then checking the quality again of these samples by FASTQC.
The FASTQC reports did not show anything out of the ordinary, so I proceeded to align the samples using the Hisat2 software.

All of these steps were done on the Linux terminal in bash.
Post alignment I went to the R software, to count reads with FeatureCounts, the command I used was as follows:

fc <- featureCounts(files = fls, 
                    annot.ext ='/home/.../Bos_taurus.ARS-UCD1.2.102.gtf' , 
                    isGTFAnnotationFile = TRUE,
                    GTF.featureType = "exon", 
                    GTF.attrType = "gene_id", 
                    useMetaFeatures = TRUE, 
                    strandSpecific = 0, 
                    isPairedEnd = TRUE, 
                    nthreads = 30)

As my samples were in two lanes, I added the read count by samples with the following command:

cont_fc2=(cont_fc[1:(ncol(cont_fc)-1)] + cont_fc[2:ncol(cont_fc)])[c(T,F)]

That done, I started to analyze the differential expression with Deseq. The results I got were few genes differentially expressed as in the example below:

out of 20328 with nonzero total read count
adjusted p-value <0.1
LFC> 0 (up): 1, 0.0049%
LFC <0 (down): 4, 0.02%
outliers [1]: 22, 0.11%
low counts [2]: 0, 0%
 (mean count <0)

Most of the analyzes were made using the “Beginner's guide to using the DESeq2 package”, by Love, M .; Anders, S. and Huber, W.
I did the test with EdgeR, but the results were similar.
To check if I did something wrong, I used the data that was used in the mentioned guide, data from the "parathyroidSE" package and also the raw RNA sequencing data mentioned in this guide so that I could do the same steps that I did with my data. I used the same parameters that I used in my data. This was done to check if something I did before the differential expression analysis, had an influence on the low results I found.
The results were as follows:

results with the ready-made data from the “parathyroidSE” package:

out of 32082 with nonzero total read count
adjusted p-value <0.1
LFC> 0 (up): 38, 0.12%
LFC <0 (down): 44, 0.14%
outliers [1]: 0, 0%
low counts [2]: 24200, 75%
(mean count <210)

results with the parameters I used in my alignment:

out of 25980 with nonzero total read count
adjusted p-value <0.1
LFC> 0 (up): 14, 0.054%
LFC <0 (down): 14, 0.054%
outliers [1]: 0, 0%
low counts [2]: 22,946, 88%
(mean count <116)

As you can see, there is a smaller amount of genes when I used the steps I did in my alignment. Can anyone help me with suggestions of what may have happened?

Thanks!



Source link