SMRTPortal环化及纠点工作

1、SMRT Portal环化及纠点工作1、 环状基因组环化目的：环状基因组环化，复制起始位点调整。1. SMRT Portal拼接基因组（HGAP_Assembly.3）,如果拼接结果很好，可以进行环化工作，则下载polished assembly fasta序列。2. Pac Bio拼接的基因组序列两端一般有重复片段，通过环化和复制起始位点调整达到去掉其中一端的重复片段的效果。使用consed工具完成环化工作。2.1 将下载的polished assembly fasta序列去掉“|quiver”。2.2 创建3个目录，这3个目录必须同级放置，进到edit_dir目录中。mkdir chrom

2、at_dir edit_dir phdball_dircd edit_dir2.3 将fasta序列转为consed识别的ace文件。fasta2Ace.perl fasta在命令行中，键入"consed"即可运行程序，程序 打开以后会弹出一个选择输入的 ace 文件的窗口：双击或者回车打开所选择的ace文件。consed主界面：2.4 检索两端是否有重复片段。点击“Search for String”，选取首端或者末尾选中2030bp，在复选框中键粘贴序列，然后点击“OK”，查找序列，如下图：如果检索到结果，界面如下：双击其中一个，会自动跳转到相应片段的相应位置。三代拼接

3、的序列首尾有overlap，那么可以进行下一步环化。2.5 确定复制起始位点的大概位置，并环化。有2种方法：第1种：根据polished assembly fasta序列，选取3K以上片段，到NCBI blast作比对，找到相应的近源菌的序列，根据基因“dnaA”的位置（该基因的product="chromosomal replication initiator protein"），获得该基因的基因序列。按照2.4步检索序列，如果有结果，如图所示：那么该基因组中复制起始位点的位置大概在 1638545 之前500bp左右。如果Search结果是“complemented”，

4、则先反转序列，再重新Search。截取片段。双击“unitig_0”，调出“Aligned Reads”窗口，点击“file”->“Export Consensus”，输入起始、终止位置，点“OK”。命名截取序列的名称。假设位于1638545 之前的片段为A，1638545 之后为B，那么调整后的顺序为 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”，如果比对区域大部分是匹配的，就连接：点击“Join contigs”。导出连接好的染色体序列，“file”->“Export Consensus”，命名结果文件。最后保存consed记录，在consed主界面，“file”->“Save assembly”

6、，退出consed。到此，基因组序列环化工作结束。第2种：如果通过NCBI blast找不到近源基因组序列，或者该近源菌的dnaA序列不太保守，在所要环化的菌中检索不到dnaA。那么只能利用GC偏移（GC skew）确定复制起始位点的大致位置。因为，在大多数细菌基因组中，前导链(leading strand)和滞后链(lagging strand)在碱基组成上存在很明显的不同前导链富含G和T，而滞后链中的A和C更多一些。打破A=T和C=G的碱基频率发生的偏移，被称之为 “AT(AT-skew)”和“GC(GC-skew)”。由于通常GC偏移比AT偏移发生的更明显，所以习惯上更多地只考虑GC偏移

7、。因为GC偏移在前导链中是正值而在滞后链中为负值，所以GC偏移值是前导链起点、终点以及转变成滞后链的信号。这使得GC偏移成为在环状染色体(circular chromosomes)中标记起点和终点的一个有用的工具。通常，当GC偏移值从负值转为正值，这一转变处位置可以认为是前导链的起点。这种方法适用于单复制起始点的细菌，对于多复制起始点或者复制起始点不明朗的菌株不太好用。先用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左右。选取一个平均数或者中位数，作为复制起始点。然后截取片段，连接，操作和第1种方法的操作相同。3. 质粒环化。如果存在质粒序列，其环化工作与上述的基因组环化相同，不过质粒可能没有dnaA基因，而且GC skew不明显，那么此时的环化就只是单纯的将首尾的重复片段除掉一份，不需要调整复制起始点。2、 基因组校点因为三代准确率和组装问题，三代拼接的序列一般存在一定的错误。使用准确率更高的illumina（二代）数据，map到

9、已环化的基因组序列上，校正基因组序列。一般使用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