多物种全基因组比对得到保守的DNA序列

查阅不同文献和教程发现获得多物种保守的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

Posted

in

by

Tags:

Comments

3 responses to “多物种全基因组比对得到保守的DNA序列”

  1. CX Avatar
    CX

    您好,您这个是用LASTZ比对吧?有详细的关于用LAST比对两物种基因组序列的吗?我看您写了一篇关于LAST安装和使用的,但是怎么建立索引呢?谢谢

    1. Wuchangsong Avatar
      Wuchangsong

      你好,last的使用可以参考http://animal.nwsuaf.edu.cn/code/index.php/Ruminantia/loadByGet?address[]=ruminantia/Documentation/documentation.php

  2. 有诺 Avatar
    有诺

    请问reflect_no_darer.txt文件是啥?

Leave a Reply to 有诺 Cancel reply

Your email address will not be published. Required fields are marked *