snakemake-workflows / rna-seq-kallisto-sleuth

A Snakemake workflow for differential expression analysis of RNA-seq data with Kallisto and Sleuth.
MIT License
66 stars 44 forks source link

fix: simpler three prime QuantSeq cutadapt setup #78

Closed dlaehnemann closed 1 year ago

dlaehnemann commented 1 year ago

I have tested this on actual QuantSeq data, and the cutadapt setup in this PR yields mostly the same results that the complicated previous three step setup recommended by Lexogen gave. It seems to be a bit more aggressive in the poly-A tail removal, but this should not affect results much. And it greatly simplifies the workflow setup.

dlaehnemann commented 1 year ago

We have now tested this on some current QuantSeq data, and the results seem to add up, even though there aren't any specific --poly-a statistics in the cutadapt log file. But using cutadapt 4.4 with Python 3.10.12 with these new recommended settings (using the standard cutadapt workflow rule):

--cores 8 -a r1adapter=AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC;min_overlap=7;max_error_rate=0.005 --minimum-length 33 --nextseq-trim=20 --poly-a

We in this instance get:

=== Summary ===

Total reads processed:              11,901,529
Reads with adapters:                 1,680,520 (14.1%)

== Read fate breakdown ==
Reads that were too short:             461,240 (3.9%)
Reads written (passing filters):    11,440,289 (96.1%)

Total basepairs processed: 1,202,054,429 bp
Quality-trimmed:              10,535,736 bp (0.9%)
Total written (filtered):  1,035,645,758 bp (86.2%)

Even when accounting for the adapter trimmed basepairs recorded in the log, this leaves us with a gap of basepairs unaccounted for of:

Total basepairs processed:  1,202,054,429 bp
Quality-trimmed:          -    10,535,736 bp
adapter trimmed:          -    48,066,270 bp
Total written (filtered): - 1,035,645,758 bp (86.2%)
============================================
--poly-a + too short          107,806,665 bp (8,97%)

As noted, these have to be basepairs removed by removing reads that are too short, plus possibly the --poly-a option. So I guess these numbers are OK and we can go ahead and merge this PR.