Genome assembly

From CSBLwiki

(Difference between revisions)

Jump to: navigation, search

Revision as of 08:12, 18 July 2010

Logbook

454 reads - de novo  +  solexa -fake reads

MIRA 사용하기
조립에 두가지 방법을 제시하고 있음
1. full de-novo 454 reads + solexa reads (총 126.9 GB 필요)
2. 454 read만으로 de-novo (2.9 GB 필요) 한 이후 solexa reads를 mapping (145.6 GB 필요)

solexa reads를 쪼개서 mapping이 가능할까?

우선 454 read만 조립
sff_extract sff_extract -l linker.fasta my_454_file.sff

fake reads -> newbler and phrap
*내가 만든 스크립트 사용
454PE-cabog -> fake reads
454SE-cabog -> 사용안함
454SE-newbler -> fake reads
454SE_PE-cabog -> 사용안함

fake reads(454PE-cabog) + fake reads(454SE-newbler) + fake reads(illu-abyss) + fake reads(illu-velvet)
1.phrap  (default) -> phrap 메모리 에러
2.newbler (-ace) -> 결과가 별로 좋지 않음, paried end 정보가 없으니 scaffold 생성도 안됨 -> 454PE reads 추가하여 scaffold 얻음, 11 -> gapRes -> 각종 에러.

*MIRA fragment로 쪼개는 스크립트 + multi contigs 적용 스크립트 만들기
잘 안됨... 
pair 정보를 넣어줘야 할텐데
만약 scaffold 파일을 쪼갤경우 n을 어떻게 처리할 것인가? 그대로 두면 엄청난 참변이...
그렇다고 그냥 contig 파일을 쪼개면 무슨 의미가 있을까?
fake reads(454PE-cabog) + fake reads(454SE-newbler) + fake reads(illu-abyss) + fake reads(illu-velvet)

 다음 step 
*cabog에 들어가는 fastq의 길이 확인 -> contig를 fake read로 만들기 -> 조립
*cabog의 contig를 fake read로 만들고 -> newbler로 조립 -> gapRes
*small assembly를 만들어서(ace 파일등) -> dupfinisher 디버깅
*phrap 으로 fake read를 조립 -> ?
*cabog 를 gapRes이 사용하도록 변경

 cabog with ace output and some options 
~/tools/wgs-6.1/Linux-amd64/bin/runCA -d SE -p SE createACE=1 unitigger=bog doToggle=1 closureOverlaps=0 closurePlacement=2 SE.frg & ~/tools/wgs-6.1/Linux-amd64/bin/runCA -d PE -p PE createACE=1 unitigger=bog doToggle=1 closureOverlaps=0 closurePlacement=2 PE.frg & ~/tools/wgs-6.1/Linux-amd64/bin/runCA -d SE_PE -p SE_PE createACE=1 unitigger=bog doToggle=1 closureOverlaps=0 closurePlacement=2 SE.frg PE.frg &

gapResolution 사용
/home/gnusnah/works/assembly_2010_7_8/gapRes/run1
~/tools/gapResolution-1_2_1/bin/runGapResolution.pl -od run1 -np 8 ../SE_PE_abyss/assembly/consed/edit_dir/454Contigs.ace.1 ../SE_PE_abyss/assembly/454Scaffolds.txt ../SE_PE_abyss/assembly/454NewblerMetrics.txt ../SE_PE_abyss/assembly/454AllContigs.fna ../SE_PE_abyss/assembly/454AllContigs.qual
~/tools/gapResolution-1_2_1/bin/stitchClosedSubProjects.pl ../../SE_PE_abyss/assembly/454Scaffolds.txt ../../SE_PE_abyss/assembly/454AllContigs.fna ../../SE_PE_abyss/assembly/454AllContigs.qual ./fakes/ ./assemInfo/gapdirs.txt my_run1
~/p-code/PModule/assembler_modules/scf2ctg.py my_run1.fasta

seqanswers에서 mira 3의 사용이 hybrid에 상당히 유효하다는 의견들이 있음
메뉴얼이 consed 못지 않게 김.

cabog 사용, read:454PE,454SE,illumina 2
만 1일째 0 단계 overlap 중, 언제 끝날지 예측 불가. cpu 사용양을 보니 190%. 몇개를 이용하는지는 알 수 없음. 0-overlaptrim-overlap 단계에서 하드디스크 용량 문제로 실패. 실패한 부분에서 무려 64GB를 차지함.

 St. Louis conversion script 제작 중 
제작 중 454 오리지널 read를 살펴보니, mate pair 정보가 들어있는 read의 경우 linker seq로 쪼갠 후 양 끝 중 어느 한쪽이 짧을 경우 정보를 버린다는 것을 알게됨.
그래서 newbler를 이용해 최소 read 길이 옵션을 조정해서 조립함. 20(default) -> 15(바꿀 수 있는 최소길이)
결과는 오히려 더 안좋아짐. 이 것은 아마도 짧은 서열은 더 많은 혼동을 주기 때문으로 생각됨
script 제작 중 qual 정보를 다루는 것이 어려워 잠시 중단

cabog 사용, read:454PE,454SE,abyss contigs
panpyro
실패 fastq를 읽는 부분은 illumina read에 맞도록 되어 있는 것으로 생각됨. 긴 read는 읽히지 않는 것 같음.

cabog 사용, read:454PE,454SE,abyss fake reads
panpyro /home/users/roh329/works/assembly_2010_7_12
실패 abyss fake reads에 알 수 없는 문제가 있음

fake qual을 만들고 fasta와 섞어서 fastq만듬
/home/gnusnah/p-code/PModule/assembler_modules/make_qual.py
/home/gnusnah/p-code/PModule/assembler_modules/make_fastq.py

cabog 사용, read:454PE,454SE,illumina
panflam
~/tools/wgs-6.1/Linux-amd64/bin/fastqToCA -insertsize 375 25 -libraryname JUN_illu -type illumina -fastq /home/gnusnah/db/genome/Eubacteria/JUN_2010_PE/s_3.1.fastq,/home/gnusnah/db/genome/Eubacteria/JUN_2010_PE/s_3.2.fastq > s_3.frg
~/tools/wgs-6.1/Linux-amd64/bin/sffToCA -libraryname PE -insertsize 3000 200 -linker titanium -output PE GE6FA8204.sff
~/tools/wgs-6.1/Linux-amd64/bin/sffToCA -libraryname SE -output SE GIST.SE.sff
~/tools/wgs-6.1/Linux-amd64/bin/runCA -d SE_PE_ILLU -p run1 unitigger=bog doToggle=1 clossurePlacement=1 PE.frg SE.frg s_3.frg

abyss contigs의 fake reads + 454 data
phrap 사용이 어려워, newbler로 조립해봄, commandline manual을 못찾아 GUI로 조립: -consed -a 50 -l 350 -ml 20
scaffold: 11->8, contigs수: 64->290, contigs총길이: 4247430->4284534

solexa reads로 만든 abyss contigs의 fake read 만들기
길이는 1.5kb, 그 이하의 contigs는 다 버려야 하나? phrap으로 조립하기 위해서는 아마도...
coverage는 얼마나? 10
/home/gnusnah/p-code/PModule/assembler_modules/make_randomread_4_illu_contig.py
45221개, 총길이 67828507의 라이브러리 만듬

phrap 사용 solexa 조립
read의 이름을 어떻게 변환? manual을 보면 "create a script which translates your read names into St. Louis", 다른 사람들이 만들어 놓은 script는 없나?

다시 addSolexaReads.perl
gnusnah@panflam:~/works/assembly_2010_7_8/SE_PE/consed/edit_dir$ addSolexaReads.perl 454Contigs.ace.1 solexa_files.fof ref.fa 
약 2시간 걸림, 또 실패
couldn't execute /home/gnusnah/tools/UW/consed/bin/consed -ace 454Contigs.ace.1 -addReads alignmentFiles100711_154311.fof -chem solexa at /home/gnusnah/tools/UW/consed/bin/addSolexaReads.perl line 170.
error_at_reading_step quality value를 읽는 과정 -> 메모리부족 -> solexa read 자체를 읽어 들이는 것은 비효율적인것으로 생각됨 -> 논문에서처럼 contigs 쪼개서 fake reads를

100711 Solexa read 변환
"." 을 N 으로 변환: cat s_3.1.fastq | perl -pi -e 's/\./N/g' > N_s_3.1.fastq

Add solexa reads to Newbler result
gnusnah@panflam:~/works/assembly_2010_7_8/SE_PE/consed/edit_dir$ addSolexaReads.perl 454Contigs.ace.1 solexa_files.fof ref.fa 
총 33분 걸림
error - 454Contigs.ace.2 file: 0 -> 하드가 100% 됐었음, 정리 후 다시 실행
다시 error - read에 포함된 "." 가 문제 - 어떻게 해결? "." 가 있는 read 삭제? 삭제할 때는 pair인 read도 함께 삭제? -> "."을 n으로 바꾸면 될지도.

run Newbler PE
runAssembly -o PE -a 50 -l 350 -g -m -ml 20 -cpu 0 -consed ~/db/genome/Eubacteria/APR_2010_PE/GE6FA8204.sff
(/home/gnusnah/works/assembly_2010_7_8/)

run Newbler SE
runAssembly -o SE -a 50 -l 350 -g -m -ml 20 -cpu 0 -consed ~/db/genome/Eubacteria/NOV_2009_SE/GIST.SE.sff
(/home/gnusnah/works/assembly_2010_7_8/)

add solexa read, doing...
under /home/gnusnah/works/assembly_2010_7_8/consed/
make dir : solexa_dir
link to fastq (2 paired end file)
make file : edit_dir/solexa_files.fof

Consed Customization
file : /home/gnusnah/.consedrc
add environment : /home/gnusnah/.bashrc

Consed Install
Consed_Install
While customizing phredPhrap, the location of polyphred should be confirmed. Polyphred is not installed. Sent request e-mail.

run Newbler SE + PE
runAssembly -o SE_PE -a 50 -l 350 -g -m -ml 20 -cpu 0 -consed ~/db/genome/Eubacteria/NOV_2009_SE/GIST.SE.sff ~/db/genome/Eubacteria/APR_2010_PE/GE6FA8204.sff
(/home/gnusnah/works/assembly_2010_7_8/)

Try Consed
gnusnah@panflam:~/works/assembly_2010_7_8/SE_PE/consed/edit_dir$ ~/tools/UW/consed/consed_linux64bit

phred
add environment : /home/gnusnah/.bashrc
PHRED_PARAMETER_FILE=/home/gnusnah/tools/UW/phred/phredpar.dat
export PHRED_PARAMETER_FILE

Softwares


Software	Version	Input	Output	Location(machine/folder)
Newbler	2.3(091027_1459)			panflam,panpyro
Phrap	0.990329(Phrap0.990329_patch)			panflam
Phrap	1.090518			panflam
Consed	090206			panflam
CABOG(celera)	6.1	sanger, 454(.sff), illumina(fastq), fastq	CABOG_output	panflam,panpyro
maq	0.7.1	ref:fasta, read:illumina, long read(not good)		panflam,panpyro
abyss [[1]]	1.2.0	454, illumina		panflam
SOAPdenovo	1.04	illumina		panflam
Corrector(soap package)	1.00	fasta,fastq		panflam
GapCloser(soap package)	1.10	fasta,fastq		panflam
MIRA		sanger,454,illumina
gapResolution		many 454 results	fasta,qual
Dupfinisher		ace file