Dinh/Dinh 2012/NOTES/2012-3-5

From ZhangLabWiki
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.