So, I am using cutadapt to remove two primers from some paired end reads.
I have 2 files for each sample, i.e. sample_R1.fastq and sample_R2.fastq.
I have primer_1 and primer_2 and both of them can be found in sample_R1.fastq and sample_R2.fastq. In other words, if I search for primer_1 and primer_2 in sample_R1.fastq, the sum of the matches is the same as the number of the reads in sample_R1.fastq. Same for R2.
To properly trim them what I do is the following, through R:
setOne <- paste("-g", primer_1, "-a", primer_1)
setTwo <- paste("-G", primer_2, "-A", primer_2)
system2(cutadapt, args = c(setOne, setTwo, "-n", 2,
"-o", outOne, "-p", outTwo,
inputOne, inputTwo, "--minimum-length", 1))
The command works fine except for one thing. On the one hand cutadapt
should remove those sequences which, after the trimming, will have a length < 1
, as specified by the command. However, some reads are removed even though they are longer than 1 and I would like to know why this happens, since I can't find a specific reference mentioning that cutadapt
removes sequences.
If I am getting this straight, cutadapt
only removes the primer sequence from the beginning of the read, then removes the equivalent characters in the quality field of the read that was trimmed and, finally, it returns the trimmed read with the quality field adjusted to the new length. If the length is lower than the threshold, the read is discarded. If it is discarded in R1, it will be discarded in R2 and viceversa.
I give an example of read that is removed, this one from R1:
@some:sample:number 1:N:0:1
ACCGTCACGTCATCAGCAGCATCGACTCCTGTTTGATCCCCGTGCCTGTACTCAAGTCAACCAGTACACGCTTTTCAAGCTGCCTTCTGCCGTATTTCACACTCAGCGCATTTGCAATTGGTGTTCTTTGCACCACATATCTATCAGTTAAGCCACCGCTACTCCGCCTACCTCGTACTGACTTAATCGACCGCCTACGCACCCTTTAAACCCAATTCGAGACCGTCAATTCAATCATACTTTAGTAAGTG
+
BBBBBBBBBBBBGGGGHEGHFHA@EGGGHHHHHHHHHGHHHHHHHHGGGGGHGHGHGGGHHHHHGDGGGFHGGGGHGGGGGGHGHGFHHHHHHHHHHGGGHHHGHGHHHHHHHGFFGGGGGGGGGGGGGGGGGFFBDFDFFFFFFFBFFFFFFFFFHGHGGHHHHHEGGGGGHHHHGHHHHHGHHHHHHHHHHHHGGGHGHHGGGGEFGGHHHHHHHHHHHHHGHHHHHHHHHHHFHFGHHHHGGGGGHHG
this one from R2:
@some:sample:number 2:N:0:1
CTGCTGCTATCACTAGCTAAAGCTTTAACGTAGGGCAATTGACGAAGGTGAAATGCATAGATATGTGATGACGGTGTAGCGTTGGCGCTGGGGTTTAAAGGGTGCGTAGTGACTTCACAGAACACCAATTGGCAGCTTACTAAAGCATGAGTACAGACGAGGTAGGCGGAATTGCGGTCGATTAAGTCAGTTGTGAAATACGGCAGCTTAACTGTTGGCAAGCGTTGTCCGGCACTGGTTTTATTCGAGGT
+
AAAAA@BFFFFFGFGCGGFDDC/CDHF2FGFHHG//<<CCC1FGHH1<1>11C<FF/0BFFFFFFF/:DEEEFFHGG0AEG/B9B./-.@9GHGEFC>EE<FHEHEGHDHGGFHGGBGGGHGGDGH0CDDGFFHHCC--.:CGCFG??GGAEGGHFHHBFFFEFB./B..;/;:/999=AAA../BB/;//////9/EFFCBDAA..990ADD=99/;/BFGGGGGGGGGFHHHHHHGHDGHHFFDHHHHG
and these are the primers that I am removing:
- ACCGTCACGTCATCAGCAGCATCGAC
- CTGCTGCTATCACTAGCTAAAGCTTTA
matching the reads I provided. Clearly, if I trim them manually, the reads length is way greater than 1 so, why are these reads removed from the final cutadapt
output file?
EDIT:
I might have solved the issue, but I asked on the cutadapt
github issues page, just to make sure.