查阅不同文献和教程发现获得多物种保守的DNA序列(编码区和非编码区)主要过程有:
1.Repeat mask:通过RepeatMasker和RepeatModeler获得
2.Pairwise alignment: 用到的软件主要有last、lastz、blastz
3.Chaining: axtChain
4.Netting: chainNet
5.Mafing:
6.Combine multiple pairwise results:
7.PhastCons: PHAST
详细步骤如下
第一步:从数据库下载重复序列屏蔽后的基因组fasta文件,对自己组装的序列可以通过Geta获得
第二步:前边介绍过last的使用,但看文献发现使用lastz的比较多,有关last和lastz的比较(last aligner is considered faster and memory efficient. It creates maf file, which can converted to psl files. Then the same following processes can be used on psl files. Different from lastz, last aligner starts with fasta files. The target genome sequence has to build the index file first, and then align with the query genome sequence.),操作上last使用起来更加简单,参数选择较少,目前还不知道两者结果的异同(服务器正在运行,结果出来更新)。lastz和blastz的不同。
lastz数据预处理:
for i in `cat 00.initial_data/reflect.txt | cut -f 1` do echo "faToTwoBit 00.initial_data/$i.genome.fasta 00.initial_data/$i.genome.2bit" done > fa2bit.list ParaFly -c fa2bit.list -CPU 19 for i in `cat 00.initial_data/reflect.txt | cut -f 1` do echo "twoBitInfo 00.initial_data/$i.genome.2bit stdout | sort -k2rn > $i.chrom.sizes" done > chrom.sizes.list ParaFly -c chrom.sizes.list -CPU 19 for i in `cat 00.initial_data/reflect.txt | cut -f 1` do mkdir ${i}PartList done for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do echo "/opt/biosoft/userApps/kent/src/hg/utils/automation/partitionSequence.pl 10000000 0 00.initial_data/$i.genome.2bit $i.chrom.sizes 1 -lstDir ${i}PartList > $i.part.list" done > query_partitionSequence.list ParaFly -c query_partitionSequence.list -CPU 18 /opt/biosoft/userApps/kent/src/hg/utils/automation/partitionSequence.pl 20000000 10000 00.initial_data/darer.genome.2bit darer.chrom.sizes 1 -lstDir darerPartList > darer.part.list grep -v PartList darer.part.list > darer.list for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do echo "grep -v PartList $i.part.list > $i.list" done > 1111.lsit ParaFly -c 1111.lsit -CPU 18 for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do echo "cat ${i}PartList/*.lst >> $i.list" done > cat_Part.list ParaFly -c cat_Part.list -CPU 18 /opt/biosoft/userApps/kent/src/hg/utils/automation/constructLiftFile.pl darer.chrom.sizes darer.list > darer.lift for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do echo "/opt/biosoft/userApps/kent/src/hg/utils/automation/constructLiftFile.pl $i.chrom.sizes $i.list > $i.lift" done > constructLiftFile.list ParaFly -c constructLiftFile.list -CPU 18 for i in `cat 00.initial_data/reflect.txt | cut -f 1` do mkdir $i for x in `cat $i.list` do y=${x/*2bit:/} echo "twoBitToFa $x $i/$y.fa" done >> twoBitToFa.list done #去除长度小于1000bp的序列 for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do for x in $i/*fa do y=${x/*-/} k=${y/.fa/} if [ $k -le 1000 ] then rm $x fi done done ParaFly -c twoBitToFa.list -CPU 80 for i in darer/*fa do for x in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do for y in $x/*fa do echo "lastz $i $y --strand=both --seed=12of19 --notransition --chain --gapped --gap=400,30 --hspthresh=3000 --gappedthresh=3000 --inner=2000 --masking=50 --ydrop=9400 --scores=/opt/biosoft/GenomeAlignmentTools/HoxD55.q --format=axt > ${x}_axt/$i.$y.axt" done >> lastz_all.list done done #lastz的运行速度太慢,若分析的物种太多没有超算不推荐使用
第三步:Chaining,将相邻的block连接起来,打分矩阵和blastz相同,gap打分改变
for i in `cat 00.initial_data/reflect_no_darer.txt | cut -f 1` do echo "axtChain -linearGap=loose -psl $i.psl darer.genome.2bit $i.genome.2bit $i.Todarer.chain" done > chain_axtChain.list ParaFly -c chain_axtChain.list -CPU 18
第四步:Netting:chainNet,对target序列确定最优比对序列。
1.首先将所有的染色体或scaffold的碱基标记未用的。
2.按打分由高到低排列,形成列表。
3.迭代:每次从列表中取出一个chain,扔掉与已经存在的chain有overlap的区域,余下的部分添加上去,如果和之前的chain有gap,标记成子集,通过这种方式形成的层级称为net。记录overlap的区域,用于下一步识别重复。
chainMergeSort $output_dir/3.chain/*.chain > $output_dir/4.prenet/all.chain chainPreNet $output_dir/4.prenet/all.chain $output_dir/$tn.sizes $output_dir/$qn.sizes $output_dir/4.prenet/all_sort.chain chainNet $output_dir/4.prenet/all_sort.chain $output_dir/$tn.sizes $output_dir/$qn.sizes $output_dir/5.net/temp.tn $output_dir/5.net/temp.qn netSyntenic $output_dir/5.net/temp.tn $output_dir/5.net/$tn.net netSyntenic $output_dir/5.net/temp.qn $output_dir/5.net/$qn.net
第五步:Mafing
netToAxt $output_dir/5.net/$tn.net $output_dir/4.prenet/all_sort.chain $output_dir/$tn.2bit $output_dir/$qn.2bit $output_dir/6.net_to_axt/all.axt axtSort $output_dir/6.net_to_axt/all.axt $output_dir/6.net_to_axt/all_sort.axt axtToMaf -tPrefix=$tn -qPrefix=$qn $output_dir/6.net_to_axt/all_sort.axt $output_dir/$tn.sizes $output_dir/$qn.sizes $output_dir/7.maf/all.maf
第六步:Combine multiple pairwise results:
roast + E=darer tree.txt ./*Final.maf all.roast.maf
Leave a Reply