SMRTPortal環(huán)化及糾點工作

1、SMRT Portal環(huán)化及糾點工作1、 環(huán)狀基因組環(huán)化目的：環(huán)狀基因組環(huán)化，復制起始位點調(diào)整。1. SMRT Portal拼接基因組（HGAP_Assembly.3）,如果拼接結(jié)果很好，可以進行環(huán)化工作，則下載polished assembly fasta序列。2. Pac Bio拼接的基因組序列兩端一般有重復片段，通過環(huán)化和復制起始位點調(diào)整達到去掉其中一端的重復片段的效果。使用consed工具完成環(huán)化工作。2.1 將下載的polished assembly fasta序列去掉“|quiver”。2.2 創(chuàng)建3個目錄，這3個目錄必須同級放置，進到edit_dir目錄中。mkdir chrom

2、at_dir edit_dir phdball_dircd edit_dir2.3 將fasta序列轉(zhuǎn)為consed識別的ace文件。fasta2Ace.perl fasta在命令行中，鍵入"consed"即可運行程序，程序 打開以后會彈出一個選擇輸入的 ace 文件的窗口：雙擊或者回車打開所選擇的ace文件。consed主界面：2.4 檢索兩端是否有重復片段。點擊“Search for String”，選取首端或者末尾選中2030bp，在復選框中鍵粘貼序列，然后點擊“OK”，查找序列，如下圖：如果檢索到結(jié)果，界面如下：雙擊其中一個，會自動跳轉(zhuǎn)到相應(yīng)片段的相應(yīng)位置。三代拼接

3、的序列首尾有overlap，那么可以進行下一步環(huán)化。2.5 確定復制起始位點的大概位置，并環(huán)化。有2種方法：第1種：根據(jù)polished assembly fasta序列，選取3K以上片段，到NCBI blast作比對，找到相應(yīng)的近源菌的序列，根據(jù)基因“dnaA”的位置（該基因的product="chromosomal replication initiator protein"），獲得該基因的基因序列。按照2.4步檢索序列，如果有結(jié)果，如圖所示：那么該基因組中復制起始位點的位置大概在 1638545 之前500bp左右。如果Search結(jié)果是“complemented”，

4、則先反轉(zhuǎn)序列，再重新Search。截取片段。雙擊“unitig_0”，調(diào)出“Aligned Reads”窗口，點擊“file”->“Export Consensus”，輸入起始、終止位置，點“OK”。命名截取序列的名稱。假設(shè)位于1638545 之前的片段為A，1638545 之后為B，那么調(diào)整后的順序為 B-A。重復做2.2，2.3步：mkdir chromat_dir edit_dir phdball_dircd edit_dircat temp1.fasta temp2.fasta > temp.fastafasta2Ace.perl temp.fasta將上述截取到的片段序列

5、cat到同一個fasta文件中，更改fasta ID，確保temp.fasta里的ID是唯一的，否則fasta2Ace.perl會報錯。通過“Search for String”檢索overlap，然后人為連接起來：點擊“Compare cont”(兩個片段的都要點擊)，將2個要連接的片段放到Align窗口中，點擊“Align”，如果比對區(qū)域大部分是匹配的，就連接：點擊“Join contigs”。導出連接好的染色體序列，“file”->“Export Consensus”，命名結(jié)果文件。最后保存consed記錄，在consed主界面，“file”->“Save assembly”

6、，退出consed。到此，基因組序列環(huán)化工作結(jié)束。第2種：如果通過NCBI blast找不到近源基因組序列，或者該近源菌的dnaA序列不太保守，在所要環(huán)化的菌中檢索不到dnaA。那么只能利用GC偏移（GC skew）確定復制起始位點的大致位置。因為，在大多數(shù)細菌基因組中，前導鏈(leading strand)和滯后鏈(lagging strand)在堿基組成上存在很明顯的不同前導鏈富含G和T，而滯后鏈中的A和C更多一些。打破A=T和C=G的堿基頻率發(fā)生的偏移，被稱之為 “AT(AT-skew)”和“GC(GC-skew)”。由于通常GC偏移比AT偏移發(fā)生的更明顯，所以習慣上更多地只考慮GC偏移

7、。因為GC偏移在前導鏈中是正值而在滯后鏈中為負值，所以GC偏移值是前導鏈起點、終點以及轉(zhuǎn)變成滯后鏈的信號。這使得GC偏移成為在環(huán)狀染色體(circular chromosomes)中標記起點和終點的一個有用的工具。通常，當GC偏移值從負值轉(zhuǎn)為正值，這一轉(zhuǎn)變處位置可以認為是前導鏈的起點。這種方法適用于單復制起始點的細菌，對于多復制起始點或者復制起始點不明朗的菌株不太好用。先用artemis繪制GC skew圖，確定復制起始點大概位置。在命令行里鍵入：art，打開Art界面。由此，可以得知復制起始點在1630K之后。python /work/xzh/TOOLS/bin/circos_tools/G

8、Ccalc.py -f unitig_0.fa -w 2000 -s 1000 > temp.gcskewmore temp.gc可以看出復制起始點在1637K1639K左右。選取一個平均數(shù)或者中位數(shù)，作為復制起始點。然后截取片段，連接，操作和第1種方法的操作相同。3. 質(zhì)粒環(huán)化。如果存在質(zhì)粒序列，其環(huán)化工作與上述的基因組環(huán)化相同，不過質(zhì)?？赡軟]有dnaA基因，而且GC skew不明顯，那么此時的環(huán)化就只是單純的將首尾的重復片段除掉一份，不需要調(diào)整復制起始點。2、 基因組校點因為三代準確率和組裝問題，三代拼接的序列一般存在一定的錯誤。使用準確率更高的illumina（二代）數(shù)據(jù)，map到

9、已環(huán)化的基因組序列上，校正基因組序列。一般使用bwa、samtools做map，得到snp、indel信息。DynamicTrim.pl -h 20 -d ./ R1.fastqDynamicTrim.pl -h 20 -d ./ R2.fastqLengthSort.pl -l 25 -d ./ R1.fastq R2.fastqsh /work/xzh/TOOLS/bin/snp_indel/run_samtools4snp_and_filterRepeat.sh R1.fastq.trimmed.paired1 R2.fastq.trimmed.paired1 ref_sequence out_prefix Tperl /work/xzh/TOOLS/bin/snp_indel/snp_caller_from_vcf.pl vcf.file > snp_indel.xlsperl /work/xzh/TOOLS/bin/script/changeSNP_mul.pl ref_sequence snp