mourisl / Rcorrector

Error correction for Illumina RNA-seq reads
GNU General Public License v3.0
63 stars 18 forks source link

Rcorrector dropping reads? (again?) #19

Open angelaparodymerino opened 4 years ago

angelaparodymerino commented 4 years ago

Hi,

I am having this issue that was reported in a closed before. I am running Rcorrector in some RNAseq samples (first two by two) and some samples (not all of them, but the majority) seem to lose reads after the run. Run it a second time on some samples (this time one by one) and still the same problem with different number of

For example: Samples before running Rcorrector:

~/data/RNAseq_PRJNA338760/FastQrawRNAseq$ grep -c "@HISEQ1" SRR4030253_1.fastq
18218121
~/data/RNAseq_PRJNA338760/FastQrawRNAseq$ grep -c "@HISEQ1" SRR4030253_2.fastq
18218121

After running Rcorrector (first time):

grep: SRR4030253_1.fastq: No such file or directory
mcv@bambi:~/data/RNAseqCorrected$ grep -c "@HISEQ1" SRR4030253_1.cor.fq
7839585
mcv@bambi:~/data/RNAseqCorrected$ grep -c "@HISEQ1" SRR4030253_2.cor.fq
7839611

After running Rcorrector (second time)

/angela/rcorrector$ perl run_rcorrector.pl -1 ~/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_1.fastq -2 ~/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_2.fastq -od RNAseqCorrected253
Put the kmers into bloom filter
/home/mcv/angela/rcorrector/jellyfish/bin/jellyfish bc -m 23 -s 100000000 -C -t 1 -o tmp_a798458599d74d3e1d510f550790024f.bc /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_1.fastq /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_2.fastq 
Count the kmers in the bloom filter
/home/mcv/angela/rcorrector/jellyfish/bin/jellyfish count -m 23 -s 100000 -C -t 1 --bc tmp_a798458599d74d3e1d510f550790024f.bc -o tmp_a798458599d74d3e1d510f550790024f.mer_counts /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_1.fastq /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_2.fastq 
Dump the kmers
/home/mcv/angela/rcorrector/jellyfish/bin/jellyfish dump -L 2 tmp_a798458599d74d3e1d510f550790024f.mer_counts > tmp_a798458599d74d3e1d510f550790024f.jf_dump
Error correction
/home/mcv/angela/rcorrector/rcorrector -od RNAseqCorrected253  -p /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_1.fastq /home/mcv/data/RNAseq_PRJNA338760/FastQrawRNAseq/SRR4030253_2.fastq -c tmp_a798458599d74d3e1d510f550790024f.jf_dump
Stored 83145603 kmers
Weak kmer threshold rate: 0.014117 (estimated from 0.950/1 of the chosen kmers)
Bad quality threshold is '#'
Processed 36436242 reads
    Corrected 41138010 bases.

~/angela/rcorrector/RNAseqCorrected253$ grep -c "@HISEQ1" SRR4030253_1.cor.fq
10861080
~/angela/rcorrector/RNAseqCorrected253$ grep -c "@HISEQ1" SRR4030253_2.cor.fq
10861067

I have no idea what is causing this. If you could help?

Thanks in advance,

'Angela

mourisl commented 4 years ago

Which version are you using? Does it drop on both read files or only one of them? Thanks.

angelaparodymerino commented 4 years ago

Thanks for the answer.

How can I know which version is it? I could not find the way to get this info.

I got Rcorrector from: https://github.com/mourisl/Rcorrector

Let me know where should I go to find the version.

Thanks again,

'Angela

mourisl commented 4 years ago

So you used "git clone" to obtain the program? You can see the version by running "perl run_rcorrector.pl -version".

angelaparodymerino commented 4 years ago

I wasn't the person installing it but I would guess he did it using git clone. I will need to ask him though.

The version I am using is Rcorrector v1.0.4

Regards,

'Angela

mourisl commented 4 years ago

It should be the newest version then. The number of missing reads is too much. Could you please show me the last few lines from "SRR4030253_*.cor.fq" respectively? And does each read file have read ids with "unfixable_error" attached? Thanks.

angelaparodymerino commented 4 years ago

Sure:

~/angela/rcorrector/RNAseqCorrected253$ tail -n 40 SRR4030253_1.cor.fq 
CCAGCTTCAAATCCTCCACAGCTTCAACGATCTTCTCCCTCGCATCTCCCCAATACAGCTTAGTCACCACAATTATCTCTTTCTGCCGTGTAGCAGTATC
+
GGGGGGGFG@GGFGGGGGFGEFGGCGGEGEDFGGFGGFGGGGGGGGEFFDGCDGGGFGEGGGCAFFFFFGFFFF?GGGGGFGGFBGFG=GAEDEFBCEE<
@HISEQ1_0063:2:44:10266:28575_forward/1 l:1804 m:2196 h:3188 cor
ATTAAAGTGATCCTCAGAATTGATTCCAGTACCTAGGCACTGGGTGCATTTCTTACAACCTTCTCCCTCACAGTCACCACAGATAAGGCTATTTACTTTA
+
GGGFGGGDBDDGDGGGGGGDGGG=D@BCB>GEGGG5F;FFGEFFDDFECAGGFFGGBF<FEFAE<FDBEF?DD95EEFEFGCFBDGGGGGGFFDE=@EDE
@HISEQ1_0063:2:44:10321:28570_forward/1 l:3854 m:4196 h:5066
AACCACTCGGGCTTCCTTAACTGCGCCGTGCTCGCTAAACAATTTCTCCAATGAAAGATTGTCAACACCCCATGCTAGGTTGCCCACGTAGACTCGGTTG
+
HGHHHHHHHHHHHHHGGHHHHHHFHHHDEFEEHBGFHHHHFHEHHHHHHEFFFEFGFFBBEEEEFHBFFF?HHCHFHHHEFHHHHHHFBFCBDCD6D1@#
@HISEQ1_0063:2:44:10432:28583_forward/1 l:20 m:43 h:62
CTTGTCTGATCCAGATAAAAGACGCCAGTATGACAATGCCGGTTTCGAGGCTGTTGAATCAGAAAGCCAGGAATTAGAACTGGATCTTTCAAGTTTAGGA
+
HFFHGHHHEH@EGGFHHHHHGEGHGHEHFEHHHHHHHHFHGF=FGGGFEDEFF<E96@>>.@>?A9.<:<C?@AB6=;B+>=@=>?ACD9<B@>@<55C<
@HISEQ1_0063:2:44:10444:28638_forward/1 l:323 m:627 h:849
CAAAGAAGGGAAGACCCTGGACTTCAACAGTACCGAGCTTTTTCAACCCAGCGGCAAATGCTCCAGCGAGACCATGAACCCTCTGTGCAATTTGTTTAAG
+
HHGGHHHHHHGHHHHHHHHHFHHHHHHHHGHHHHHHHHHHHHHHHHHHHHHHHHHFHFHHHHFHHHHFDHBHHHFHHHHHGHHHGFGHHECFHHFEBGHH
@HISEQ1_0063:2:44:10301:28652_forward/1 l:1 m:26 h:108
GTGTAAAGGGGGGTTGAATGTGGATCTGGATTTGTAGATGGAATTTTCCCGGTTTTGAGTCCAGTTGCTCTAGACTTTCTTGACTTTGTACAGAAAGAAT
+
HHHHHHHHHHHHHHHGFGGFEDGGGHHHHFHHHGHGHEHHHHHHAHHHHHFFEHHFEFEFHHEHHFHHHHHHBEHHHHHFH@HHHGGFECFHFFFHFHH;
@HISEQ1_0063:2:44:10482:28651_forward/1 l:5 m:7 h:13 cor
GTCATGATCATCATTCAATAGTTTGGCAGATCGCCGAGGATTACTCTCATTGTGAATTTCTTATCCGTCTCCCAGGATTCGGCCCAATGCCCGCATTCCG
+
EHHGHHGHHHHHHH.GGGGFGHHHFHHFHHHHHGHHFHHFBHFHHHHHHHHFFHFHHHHHHHHFHHFBHFECFFCGCEFEHEDHHBEC9=FCCCA?DCCA
@HISEQ1_0063:2:44:10421:28691_forward/1 l:37514 m:77650 h:126177
ACGGCTCCCATCAAGATAACTTGGCAAGCCCAGATTGCCAATATGCTTTGCGCATGTACCAAGCTTGGGTTCCCCAAGTAGTCCAACCCACCTTCGCTGA
+
EEDDEEDEEBDD=DCDDDC:@@>@58@7@@EE@@EECEEBEDEEEDDEDCD:A:B<6>BB9++34;7><.=<?@@9DBDDD4DDBDBDD4D?7DACD<EB
@HISEQ1_0063:2:44:10403:28696_forward/1 l:6529 m:8230 h:9709
CTTCCGTCAATTCCTTTAAGTTTCAGCCTTGCGACCATACTCCCCCCGGAACCCAAAGACTTTGATTTCTCATAAGGTGCCGGCGGAGTCCTAAAAGCAA
+
GHHEHFHFGGDHHGHB?EEEHHHFGGBDGFGHGHGGHHEHGG?GGGGEGDEGGGFGHGEHG=GG>G@FGFGFGFBEBAE8ECDEEAAD(AD@=AD?=0?A
@HISEQ1_0063:2:44:10325:28708_forward/1 l:474 m:566 h:673 cor
ATCGTTGTAAGCTTGGCGATGGAGGTGTTCGGATTCGAGAATTACGCCATCGCAGTCGAATATGAGGGCTTGAAGGGCCCGGGAAGAAGCGGAGACGGTG
+
@EEEE@EEEEA5A7A4><><CDCDBCD=CDFBF8DFEFE8CCD8BEDBEEDCC,D;?:A;4;22:>DDCDE@EE?DCDDD@A;?BD-DB@=@########
@HISEQ1_0063:2:44:10473:28715_forwa

and

~/angela/rcorrector/RNAseqCorrected253$ tail -n 40 SRR4030253_2.cor.fq 
@HISEQ1_0063:2:44:9846:28722_reverse/2 l:52 m:116 h:130
CTTCTTTCGGATACAATAAACACATTTGTTATTTTTAATGTCAACCACAAATAAGTCTGACACCTAAATCTTGGGCCTCAGTGAACTACTAAAAGGGTGT
+
HHHHHHHHHHHHGHGHHHHHHHGHHHHHHHHHHHHHEGHGHHHHHHHHHHGGGGGHHGHHHHHHHHHHGHHHHEHHHHHHHHHFHHHHHHHHBGHHFDHF
@HISEQ1_0063:2:44:10080:28507_reverse/2 l:1116 m:1427 h:1879
TCACGAGCAATTTGTGTTACGAGAAATTTGTTTACCTGCAAGTTTGATGATGGCTCTAGTGTCTGGAGGAAGGTCGTCGTCGCTCAATCCGAATGCTCCG
+
HHHHHHHHHHHHHHHHHHGHHFHHGHHHHHHHHHHFHHHHHHDHHHHHHHHHHHHHGHHHHGHHHGCFGAFFGAFFCCBBDD@BCEAAFEDDBCE?DDFD
@HISEQ1_0063:2:44:10127:28510_reverse/2 l:49 m:125 h:213 cor
TGGAAACCTACCTGCTGAAGTGGCCAAAGACATTCATTGTTGTTTCTCATGCCCGAGAGTTCTTAAATACAGTTGTCACAGACATTCTCCATCTCCAAGG
+
@6?<ABCCDA=CA=D@:?A@A6@B@>>744>9>>==DDDDFEFGGGHGGGGDFE=:<@<4CDCC=FBDGFEFF@B;<@C4GF@GFE5BE###########
@HISEQ1_0063:2:44:10075:28557_reverse/2 l:917 m:1929 h:2309
CACAGGAGGCTGTGGTTAATGCCAGACCCTTCTCCGAAACGCTAGTCGAAGTTCTCTACAACATAAATGAGCAGCTCCAGACAGAAGATGTCGATATTCC
+
IGGIGF;DEEBEEEEEEDGFF9EFF?;:@=EEG4ED<BEE?:DE?>ECDDE?EEE@A=@@FBD@ECD;BDBFGG.EAA9DDDA#################
@HISEQ1_0063:2:44:10047:28575_reverse/2 l:256 m:307 h:390
CTTCCACCCTCAAAAATGCCCTTCTTTCTTCTCTCTTCCTCCAGTAGCCAGCCTTAGATCCCCCAACTTTTACATTGCTTCCGCTCTTCGTACTGGTTCC
+
HHHHHEHHHHHHDHBGFHCEHBHHHHHHHHHHFHHHHHHHHHBFHFE=IDEEAEE;@BCAGGG@;;ECDEGGFFGHHHDFBEF?ECFFGB:>@7@>6>C@
@HISEQ1_0063:2:44:10058:28632_reverse/2 l:2688 m:3681 h:4776 cor
GGCAGTTTTCAGGTCCGCTTGACAGGATCAAAGTTAGAGACACCAATAACCACGATGATAATGGCGGAGACCACAATACCGCCAGCACCGATGCTAAGCA
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHHHHHHHHHGHHFHHHHHHHHHHEHHHGHHHHHHHAE2BDABDB7C7A####################
@HISEQ1_0063:2:44:10116:28636_reverse/2 l:6 m:27 h:31
CCAGAGGATTTCTTCAAAATCACCAAGGAAAGAGCCTGCCTAGTGAAGTGGGCACCACAGAAGAAGGTTCTGGCACATCCATCAGTGGGAGGATTTTTCA
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHEFHHHHEHHEHHEHHHHHHHHHHHHEHHGHHEFGHFHHHHHFHEHHFHHE.EEDFEEA
@HISEQ1_0063:2:44:10196:28636_reverse/2 l:6645 m:8606 h:10217
GCCACAACCAGAACTGGGAAGTGAGAGCAGGAGGAGTGGGGAAGTGAGCCAGCTCAGTTCCGATGCTCGAGAAAAACAGCCCGGACAAGCTGTTCCCGTC
+
8CCB=8B<CD@D<DAFCECDD;C88><<<@EEEC*C3B@AEEEEDFD?6B29495;?=>=B@=@DBA@;>EEEC4:39<5EC;7?E8EE.>C?=AE4CE:
@HISEQ1_0063:2:44:10219:28703_reverse/2 l:2454 m:4650 h:5585
CTCTGTTCCATTGCTACCAAAAATAATGGAATTAGTTGGACTCGGATACACAGGATGGTTTGTGTACCGTTACCTCCTTTTCAAGTCGAGTAGAAAAGAA
+
HHHHHHHHHHHHHHHHHHHHGHHHHGHHHHHHHHHHHHHHHHHHHHHHHHHHHHEHHH<HD@DDA<>A>:D<@AA7A5ACD99A9@A5=>ABC<CDE3?B
@HISEQ1_0063:2:44:10164:28710_reverse/2 l:120 m:495 h:633
ATTCACTCAATTTCTCTCTCTAAAAATCAAACACAATAAACCCCTTCAACTCCTCCGCAAGTACAATGGCAGTGTCTCAGACCAAAGCTGTGCTCGTGAC
+
FHHHHHH

Definitely they look uncomplete.

Let me know what you think.

thanks in advance!

'Angela

mourisl commented 4 years ago

The truncation seems not an issue from Rcorrector. It feels more like the output is terminated by the operating system. Could you please check what is the remaining disk quota on your system?