Dinh/Dinh 2012/NOTES/2012-3-5
Jump to navigation
Jump to search
GM12878[edit]
- Group reads into chromosomes by initial mapping before calling ASM with Robert's pipeline.
- First, I tried to simply map using less and bisReadMapperSE_outputSOAP.pl
less s_*_sequence.txt.gz | ./bisReadMapperSE_outputSOAP.pl GM12878 101 > map.log&
- However, SOAP outputs an error at about ~40 million reads and stops.
- I deleted the log file, but one type of error is incorrect sequence size, and it printed the read name. The other one is unexpected end of file.
- Noi mentioned that with cat she could not create a file of larger than ~20 GB. So I think the same issue is with less. I was suspicious of the incorrect sequence size error and the unexpected EOF messages as well, so I decided to split the reads into 10 million read files with gzip -dc and split:
nohup gzip -dc s_7_1_sequence.txt.gz | split -l 40000000 - s_7_1. & nohup gzip -dc s_7_2_sequence.txt.gz | split -l 40000000 - s_7_2. & nohup gzip -dc s_8_1_sequence.txt.gz | split -l 40000000 - s_8_1. & nohup gzip -dc s_8_2_sequence.txt.gz | split -l 40000000 - s_8_2. &
- When checking the read which SOAP indicated as having incorrect read size, everything looked normal for that read.
- The reads were mapped as single-end with mapping rate of ~78%, but since this is whole genome, mapping as PE would have been a little better.
- GO.map.sh
for f in s_7_1.a* do ./bisReadMapperSE_outputSOAP.pl GM12878_$f 101 < $f >>map.$f.LOG done; for f in s_7_2.a* do ./bisReadMapperSE_outputSOAP.pl GM12878_$f 101 < $f >>map.$f.LOG done
- GO.2map.sh
for f in s_8_1.a* do ./bisReadMapperSE_outputSOAP.pl GM12878_$f 101 < $f >>map.$f.LOG done; for f in s_8_2.a* do ./bisReadMapperSE_outputSOAP.pl GM12878_$f 101 < $f >>map.$f.LOG done
- All the map*.LOG files were checked to make sure SOAP had no error message. The range of mapping rates were between 10-39%. This appeared odd, but perhaps some cluster of reads were coming from a poorly imaged part of the flowcell (bubble?)
- All the mapping were ok except for the LAST reads sets:
Error 1: Length Error: incompitable seq and qual length HWI-ST1001:8:2307:15776:79017#0/3|TTTTTTTTTTTTAAAAAAAAATTAAAATTTAAAATTTTTTTTAAGTTTTTTTTTTTTTTTTTTTTAGTTTATTAGAGTTGAATGTTTATAATAGGTTTAGG Error 2: File Error: unexpected feof
- I went ahead with the next steps.
- For s_7_1, s_7_2, s_8_1, and s_8_2 in parallel, the reads were combined and sorted, then filtered (clean up double hits), and then split into chromosomes:
for f in aa ab ac ad ae af ag ah ai aj ak al am an ao ap aq ar as at au av do sed 's/chr/REV:chr/g' GM12878_s_7_1.$f.rev.soap.out >> GM12878_s_7_1.$f.fwd.soap.out sort -k1,1 < GM12878_s_7_1.$f.fwd.soap.out > combined.sorted.GM12878_s_7_1.$f.soap.out ./filterHits.pl < combined.sorted.GM12878_s_7_1.$f.soap.out | ./soap2splitReads.pl GM12878.s_7_1 done
- Finally, the 4 separate files from s_7_1, s_7_2, s_8_1, s_8_2 were combined:
for f in chr1 chr2 chr3_21 chr4_22 chr5_19 chr6_20 chr7_18 chr8_17 chr9_16 chr10_15 chr11_14 chr12_13 chrX_Y do cat GM12878.s_7_2.$f.fastq GM12878.s_8_1.$f.fastq GM12878.s_8_2.$f.fastq >> GM12878.s_7_1.$f.fastq rm GM12878.s_7_2.$f.fastq GM12878.s_8_1.$f.fastq GM12878.s_8_2.$f.fastq done
- Get the last reads from fastq file:
cat s_7_1_2306_qseq.txt s_7_1_2307_qseq.txt s_7_1_2308_qseq.txt | /home/kunzhang/bin/qseq2fastq.pl > /media/Ext12T/DD_Ext12T/HL_WGBS_map/Sequences/last_s_7_1.av & cat s_7_2_2306_qseq.txt s_7_2_2307_qseq.txt s_7_2_2308_qseq.txt | /home/kunzhang/bin/qseq2fastq.pl > /media/Ext12T/DD_Ext12T/HL_WGBS_map/Sequences/last_s_7_2.av & cat s_8_2_2306_qseq.txt s_8_2_2307_qseq.txt s_8_2_2308_qseq.txt | /home/kunzhang/bin/qseq2fastq.pl > /media/Ext12T/DD_Ext12T/HL_WGBS_map/Sequences/last_s_8_2.av & cat s_8_1_2306_qseq.txt s_8_1_2307_qseq.txt s_8_1_2308_qseq.txt | /home/kunzhang/bin/qseq2fastq.pl > /media/Ext12T/DD_Ext12T/HL_WGBS_map/Sequences/last_s_8_1.av &
- Combine the reads
2023 ./combine2files.pl s_7_1.av last_s_7_1.av HWI-ST1001:7:2307:11214:103165 > correct.s_7_1.av & 2026 ./combine2files.pl s_7_2.av last_s_7_2.av HWI-ST1001:7:2307:15372:125378 > correct.s_7_2.av & 2027 ./combine2files.pl s_8_1.av last_s_8_1.av HWI-ST1001:8:2306:5020:194963 > correct.s_8_1.av & 2028 ./combine2files.pl s_8_2.av last_s_8_2.av HWI-ST1001:8:2307:15776:79017 > correct.s_8_2.av &
- How to combine the files:
#!/usr/bin/perl -w use strict; my $f1 = $ARGV[0]; my $f2 = $ARGV[1]; my $rname = $ARGV[2]; open(INFILE, "$f1")||die("Error opening file\n"); while(my $line1 = <INFILE>){ my $line2 = <INFILE>; my $line3 = <INFILE>; my $line4 = <INFILE>; last if($line1 =~ m/$rname/); print $line1,$line2,$line3,$line4; } close(INFILE); my $on = 0; open(INFILE, "$f2")||die("Error opening file\n"); while(my $line1 = <INFILE>){ my $line2 = <INFILE>; my $line3 = <INFILE>; my $line4 = <INFILE>; $on = 1 if($line1 =~ m/$rname/); print $line1,$line2,$line3,$line4 if($on == 1); } close(INFILE);
- File have the correct number of lines:
60864000 correct.s_7_1.av 60864000 correct.s_7_2.av 55590316 correct.s_8_1.av 55590316 correct.s_8_2.av
Conclusion[edit]
- Initially, we had ~10millionx22x4 ~ 800 million reads (200 million per lane). After mapping and cleaning, there were ~12x13x4 ~ 600 million uniquely mappable reads.
- For the next WGBS, I would try to map as PE
- I think splitting the reads was a good idea as we can make combine, sort, and filtering/cleaning step go in parallel. I could have parallelize this step even more.