Editing
Dinh/Dinh 2012/NOTES/2012-3-5
Jump to navigation
Jump to search
Warning:
You are not logged in. Your IP address will be publicly visible if you make any edits. If you
log in
or
create an account
, your edits will be attributed to your username, along with other benefits.
Anti-spam check. Do
not
fill this in!
==GM12878== * 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== * 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.
Summary:
Please note that all contributions to ZhangLabWiki may be edited, altered, or removed by other contributors. If you do not want your writing to be edited mercilessly, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource (see
ZhangLabWiki:Copyrights
for details).
Do not submit copyrighted work without permission!
Cancel
Editing help
(opens in new window)
Navigation menu
Personal tools
Not logged in
Talk
Contributions
Create account
Log in
Namespaces
Page
Discussion
English
Views
Read
Edit
View history
More
Search
Navigation
Main Page
Current events
Recent changes
Random page
Investigators
Matt Cai
Song Chen
Eric Chu
Dinh Diep
Elizabeth Duong
Shicheng Guo
Alan Fung
Daniel Jacobsen
Blue Lake
Huy Lam
Alice Li
Andrew Richards
Brandon Sos
Chris Wei
Yan Wu
Kun Zhang
Tools
What links here
Related changes
Special pages
Page information