0.代碼來源
文章參考:
《Large-scale ruminant genome sequencing provides insights into their evolution and distinct traits》
DOI: 10.1126/science.aav6202
后面簡稱'原文'
本文很貼心的給出了整個保守元件(CE, Conserved Element)的鑒定流程以及所用到的腳本.
https://zenodo.org/record/2549147
如果用到這個代碼，請注明代碼來源并引用這篇文章
- 軟件的安裝
參見: http://compgen.cshl.edu/phast
2.多基因組比對
這里也可以參考以下兩個鏈接
https://biochen.org/cn/blog/2021/%E5%A2%9E%E5%BC%BA%E7%89%88%E4%B8%A4%E5%9F%BA%E5%9B%A0%E7%BB%84%E6%AF%94%E5%AF%B9/
https://biochen.org/cn/blog/2021/%E5%A4%9A%E5%9F%BA%E5%9B%A0%E7%BB%84%E6%AF%94%E5%AF%B9/
假設已經(jīng)拿到了最后的maf文件
這里只演示一條染色體的結(jié)果
3.1 mod文件的獲取
PS: 原文獲取mod的做法較為繁瑣
具體見本文末尾 4. 一節(jié), 有完整的復現(xiàn)
這里先展示一種較為簡單的方法,
此處參考: https://biochen.org/cn/blog/2021/%E4%BD%BF%E7%94%A8PHAST%E5%88%86%E6%9E%90%E4%BF%9D%E5%AE%88%E6%80%A7/
以及官方給出的pipeline (真的比原文的做法簡單太多):
http://compgen.cshl.edu/phast/phyloP-tutorial.php
phyloFit -i MAF Chr1.maf --tree "(((((A, B), (C, D)), E), F), G)"
#輸出文件phyloFit.mod
#用作后續(xù)分析
輸出文件phyloFit.mod以用作后續(xù)分析
3.2 CE的鑒定
phastCons Chr1.maf phyloFit.mod --target-coverage 0.25 --expected-length 12 --rho 0.4 --msa-format MAF --seqname Chr1 --most-conserved Chr1.CE.bed > Chr1.conserved.score.wig
輸出文件1是Chr1.CE.bed, 這個就是保守元件的bed文件, 還有一個Chr1.conserved.score.wig是CE每個bp的保守性得分
其實做到這里就已經(jīng)算是完成了, 但是可以再多做一點
比如我的進化樹包含了 A - G總共7個物種，但是A-D才是我的研究物種, 我想把CE過濾到A-D上
按照原文的做法, 有兩種過濾方式, 首先把長度小于20bp的CE都過濾掉(這一步驟本文就略去)
第一是CE的 EFG三個外群物種上是沒有序列可以比對到的
或者說只有A-D也就是研究物種中才能找到的，研究物種特有的CE
第二是CE在外群物種中也有, 但是在研究物種中是 '保守的' 或者 '加速進化的' (基于看不懂的LRT test一大堆, 參考http://compgen.cshl.edu/phast/help-pages/phyloP.txt)
3.3 過濾CE (具體用到的腳本在本文末尾或者從zendo上自己下載)
3.3.1 首先是得到第一種CE:
修改腳本maf_bed_filter_Type1.pl
把這個腳本的標注出來的這一行, 修改為maf文件中外群物種的具體名稱
修改我高亮的地方
perl maf_bed_filter_Type1.pl Chr1.CE.bed Chr1.maf ./
#會生成work.out， 這個就是外群沒有，研究物種特有的CE，但是格式不是bed格式
#下面就是把work.out轉(zhuǎn)為bed文件
perl work_out2bed.pl work.out Chr1 ./
#輸出文件是Chr1.bed, 這就是第一種CE了
#做一下長度過濾, 低于20bp的扔掉
for i in {1..28}
do
awk '{if($3-$2>19)print}' Chr${i}.bed > Chr${i}.tmp && rm Chr${i}.bed && rename .tmp .Type1.bed *
done
3.3.2 其次是第二種CE
修改腳本01.convertMaf2List.pl
第一個紅色箭頭修改為maf文件作為參考基因組的物種的具體名稱
第二個紅色箭頭修改為maf文件的所有物種的名稱
修改紅色箭頭所示地方
修改mod文件 phyloFit.mod
在研究物種的祖先節(jié)點處加上標識, 這里表示我把我的研究物種的祖先節(jié)點命名為'Node'
修改紅色箭頭所示的地方
文件修改完畢，開始跑后續(xù)
perl 01.convertMaf2List.pl Chr1.maf Chr1
#輸出的Chr1.maf.lst用作后續(xù)
phyloP --subtree Node --msa-format MAF --features Chr1.CE.bed -method LRT --mode CONACC phyloFit.mod Chr1.maf > Chr1.phyloP.bed
#--subtree Node 這里就是前面提及的研究物種的祖先節(jié)點
#--features Chr1.CE.bed 最開始鑒定出來的CE
#-method LRT 大概是用了什么統(tǒng)計檢驗, 但我看不懂
#--mode CONACC 這個模式下除了能看研究物種相比于外群物種保守, 還能看研究物種相比于外群物種快速進化的
#最后輸出Chr1.phyloP.bed
輸出文件Chr1.phyloP.bed長這樣:
Chr1.phyloP.bed
最后一列如果小于0, 表示研究物種相比于外群物種快速進化的CE
最后一列如果大于0，研究物種相比于外群物種保守的CE
最后再卡一下顯著性閾值和做一下p值校正, 這里就不再贅述
head -n 1 Chr10.phyloP.bed > head
for i in {1..28}
do
awk '{if($9<0.05 && $9>0)print}' Chr${i}.phyloP.bed > Chr${i}.tmp && cat head Chr${i}.tmp > Chr${i}.phyloP.conservation.bed && rm Chr${i}.tmp
done
for i in {1..28}
do
awk '{if($9>-0.05 && $9<0)print}' Chr${i}.phyloP.bed > Chr${i}.tmp && cat head Chr${i}.tmp > Chr${i}.phyloP.acceleration.bed && rm Chr${i}.tmp
done
for i in {1..28}
do
awk '{if($3-$2>19)print}' Chr${i}.phyloP.acceleration.bed > Chr${i}.tmp && cat head Chr${i}.tmp > ../Acceleration/Chr${i}.acceleration.CE.bed && rm Chr${i}.tmp
done
for i in {1..28}
do
awk '{if($3-$2>19)print}' Chr${i}.phyloP.conservation.bed > Chr${i}.tmp && cat head Chr${i}.tmp > ../Type2/Chr${i}.Type2.bed && rm Chr${i}.tmp
done
#####
#使用wigToBigWig將wig文件轉(zhuǎn)換為二進制的bigWig文件后，就能在Genome Browser中展示了
#wigToBigWig in.wig chrom.sizes out.bw
for i in {1..28}
do
/data/01/user164/software/wigToBigWig Chr${i}.CE.wig Golani.Chr.v4.length Chr${i}.CE.bw
done
#假設給定一條序列比如lncRNA，我們想計算它的平均phastCons分數(shù)
#可以使用bigWigAverageOverBed實現(xiàn)
#bigWigAverageOverBed in.bw in.bed out.tab
for i in {1..28}
do
/data/01/user164/software/bigWigAverageOverBed Chr${i}.CE.bw /data/01/user164/workspace_for_whole_object/spalax_SV/10.CE/Ruminant.Pipe.9spp/CE.bed/Chr${i}.CE.bed Chr${i}.tmp;
awk '{print $1"\t"$2"\t"$4"\t"$5}' Chr${i}.tmp > Chr${i}.CE.score; rm Chr${i}.tmp;
done
cat *.CE.score | sort -Vk 1 > all.CE.score
awk '{if($2>19)print}' all.CE.score > all.CE.dayu20bp.score
4.原文的做法
4.1 4d.nonconserved.mod文件的獲取
其實和我前文中相比沒有太大改變, 只不過原文中用了4d位點
那么首先從maf文件中提取4D位點:
這里參考了: http://compgen.cshl.edu/phast/phastCons-HOWTO.html
#依然是使用phast包內(nèi)置的軟件
msa_view Chr1.maf --in-format MAF --4d --features Chr1.gff > Chr1.4d-codons.ss
msa_view Chr1.4d-codons.ss --in-format SS --out-format SS --tuple-size 1 > Chr1.4d-sites.ss
#對每條染色體重復這個操作
for i in {1..10}
do
msa_view split_maf/Chr/Chr${i}.maf --in-format MAF --4d --features gff/Chr${i}.gff > 4d.site/Chr${i}.4d.ss
msa_view 4d.site/Chr${i}.4d-codons.ss --in-format SS --out-format SS --tuple-size 1 > 4d.site/Chr${i}.4d-sites.ss
done
msa_view --unordered-ss --out-format SS --aggregate go,ca,ju,ga,Mus,Rat,Hamster,Squirrel,Rabbit Chr1.4d-sites.ss,Chr2.4d-sites.ss,Chr3.4d-sites.ss,Chr4.4d-sites.ss,Chr5.4d-sites.ss,Chr6.4d-sites.ss,Chr7.4d-sites.ss,Chr8.4d-sites.ss,Chr9.4d-sites.ss,Chr10.4d-sites.ss > all-4d.sites.ss
phyloFit --tree species.tree -i SS all-4d.sites.ss --out-root 4d.nonconserved
#用上述流程即可
#以上是標準的流程, 由于我的maf文件存在一些問題
#故進行修改
for i in {1..28}
do
msa_view split_maf/Chr/Chr${i}.maf --in-format MAF --4d --features gff/Chr${i}.gff > 9spp.4d.site/Chr${i}.tmp1
msa_view 9spp.4d.site/Chr${i}.tmp1 --in-format SS --out-format SS --tuple-size 1 > 9spp.4d.site/Chr${i}.tmp2
msa_view 9spp.4d.site/Chr${i}.tmp2 --in-format SS --out-format FASTA > 9spp.4d.site/Chr${i}.tmp3
samtools faidx 9spp.4d.site/Chr${i}.tmp3
for c in {go,ca,ga,ju,Mus,Rat,Hamster,Squirrel,Rabbit}
do
samtools faidx 9spp.4d.site/Chr${i}.tmp3 $c >> 9spp.4d.site/Chr${i}.tmp4
done
mv 9spp.4d.site/Chr${i}.tmp4 9spp.4d.site/Chr${i}.4d-sites.fasta
msa_view 9spp.4d.site/Chr${i}.4d-sites.fasta --in-format FASTA --out-format SS > 9spp.4d.site/Chr${i}.4d-sites.ss
done
msa_view --u
nordered-ss --out-format SS --aggregate go,ca,ju,ga,Mus,Rat,Hamster,Squirrel,Rabbit Chr1.4d-sites.ss,Chr2.4d-sites.ss,Chr3.4d-sites.ss,Chr4.4d-sites.ss,Chr5.4d-sites.ss,Chr6.4d-sites.ss,Chr7.4d-sites.ss,Chr8.4d-sites.ss,Chr9.4d-sites.ss,Chr10.4d-sites.ss > all-4d.sites.ss
phyloFit --tree species.tree -i SS all-4d.sites.ss --out-root 4d.nonconserved
##
這一套操作下來, 得到文件: 4d.nonconserved.mod
4.2 estimating rho
for i in {1..10}
do
mkdir Chr${i} && cd Chr${i}
msa_split Chr${i}.maf --in-format MAF --windows 1000000,0 --out-root split --out-format SS --min-informative 1000 --between-blocks 5000 #切割maf文件
done
#以單條染色體為例, 每條染色體重復以下操作
ls *.ss > ss.list
for ss in `cat ss.list`;do phastCons --estimate-rho $ss --no-post-probs $ss 4d.nonconserved.mod ;done
#4d.nonconserved.mod就是4.1的輸出文件
ls *.cons.mod >all.cons.mod.list
ls *.noncons.mod > all.noncons.mod.list
phyloBoot --read-mods all.cons.mod.list --output-average Chr1.ave.cons.mod
phyloBoot --read-mods all.noncons.mod.list --output-average Chr1.ave.noncons.mod
###
20250328更新
phyloBoot有些版本好像不能識別mod.list文件
所以得手動把所有的mod文件都放進命令行中：
phyloBoot --read-mods mod1,mod2,mod3... --output-average Chr1.ave.cons.mod
phastCons --most-conserved Chr1.CE.bed --score Chr1.maf Chr1.ave.cons.mod,Chr1.ave.noncons.mod \> Chr1.wig
#以單條染色體為例, 每條染色體重復以上操作
#如此原文的CE鑒定部分就結(jié)束, 后面對CE進行分類/過濾與本文的 3.3 一節(jié)相同
#區(qū)別在于用到的mod文件, 原文是用的4d位點生成的4D.nonconserced.mod
phyloP --subtree Node --msa-format MAF --features 1.HCE.bed --method LRT --mode CONACC 4D.nonconserced.mod 1.maf > Chr1.phyloP.bed
區(qū)別在于用到的mod文件, 原文是用的4d位點生成的4D.nonconserced.mod
phyloP --subtree Node --msa-format MAF --features 1.HCE.bed --method LRT --mode CONACC 4D.nonconserced.mod 1.maf > Chr1.phyloP.bed
5.原文中提供的腳本
vim maf_bed_filter_Type1.pl
use warnings;
use strict;
my $bed_file = $ARGV[0];#"/public/home/linzeshan/project/01.conserve/Type_one_re_fileter/new_method/work/29.bed";
my $maf_file = $ARGV[1];#"../data/maf2list/29/29.maf.lst";
my $work_dir = $ARGV[2];
my @arr = ("Mus","Rat","Hamster","Squirrel","Rabbit"); #外群物種
my $species = {};
foreach my $arr(@arr){
$species->{$arr} = 1;
}
open BED,$bed_file or die $!;
open MAF,$maf_file or die $!;
my $word = <MAF>;
my @first_line = split(/\t/,$word);
my $num = @first_line;
my @locate_arr;
for(my $i = 0;$i<$num;$i++){
push (@locate_arr,$i) if (exists ($species->{$first_line[$i]}));
}
foreach my $j(@locate_arr){
print $j."\t";
}
open OUT,">$work_dir/work.out" or die $!;
open ERR,">$work_dir/work.err" or die $!;
while(<BED>){
my @array = split(/\t/,$_);
my $pause = "nn";
for(my $i=$array[1];$i<=$array[2];$i++){
next if $i == 0;
my $chu = 1;
while(1){
my $nt;
if ($pause ne "nn"){
$nt = $pause;
$pause = "nn";
}
else {
$nt = <MAF>;
}
my @nt = split(/\t/,$nt);
if ($i>$nt[1]){
next;
}
elsif($i == $nt[1]){
foreach my $pocky(@locate_arr){
if ($nt[$pocky] ne "-"){
$chu = 0;
}
}
last;
}
else{
$pause = $nt;
print ERR "worng bed $i maf $nt[1]\n";
last;
}
}
print OUT "$i\t+\n" if $chu == 1;
print OUT "$i\t-\n" if $chu == 0;
}
}
vim work_out2bed.pl
my $file = $ARGV[0];
my $chr = $ARGV[1];
my $work_dir = $ARGV[2];
open IN,$file or die $!;
open OUT,">$work_dir/$chr.bed" or die $!;
my $over = 1;my $start=0;my $end=0;my $s_e = 1;
while(<IN>){
chomp;
my @array = split(/\t/,$_);
#print $array[0]."\t"."..".$array[1]."..\n";last;}=p
if ($array[1] eq "-"){
$end = $start if $end == 0;
print OUT "$chr\t$start\t$end\n" unless $start == 0;
$start = 0;$end = 0;
$over = 1;
next;
}
else{
if ($over == 1){
$start = $array[0];
$end = $array[0];
$over = 0;
}
else{
if ($array[0] > $end+1){
$end = $start if $end == 0;
print OUT "$chr\t$start\t$end\n" unless $start == 0;
$start = 0;$end = 0;
$start = $array[0];
$end = $array[0];
}
else{
$end = $array[0];
}
}
}
}
vim 01.convertMaf2List.pl
#! /usr/env perl
use strict;
use warnings;
my $maf=$ARGV[0];
my $output=$ARGV[1].".maf.lst";
my $ref="golani";
my @species=("go","ga","ca","ju","Mus","Rat","Hamster","Jac","Squirrel");
open O,">$output";
open E,">no_ref.maf";
#open O,"| gzip - > $output";
my @head=("chr","pos",@species);
my $head=join "\t",@head;
print O "$head\n";
my $control=0;
my $content="";
open I,"cat $maf |";
while(<I>){
next if(/^#/);
if(/^a\s+score=/){
$content=$_;
while(<I>){
if(/^s\s+/){
$content.=$_;
}
else{
last;
}
}
&analysis($content);
}
# last if($control++>=0);
}
close I;
close O;
sub analysis{
my $content=shift;
my @line=split(/\n/,$content);
my $head=shift @line;
$head=~/score=([\d\.]+)/;
my $score=$1;my $all_strand;
# print "$score\n";
my %pos;
my $isRefFound=0;
my $ref_chr="NA";
foreach my $line(@line){
my @a=split(/\s+/,$line);
my ($s,$chr,$start,$alignment_length,$strand,$sequence_length,$sequence)=@a;my $all_strand = $strand;
$chr=~/^([^\.]+)\.(.*)/;
my $species=$1;
$chr=$2;
if($species eq $ref){
$ref_chr=$chr;
$isRefFound=1;
my @base=split(//,$sequence);
if($strand eq "+"){
my $pos=$start;
for(my $i=0;$i<@base;$i++){
if($base[$i] ne "-"){
$pos++;
$pos{$i}=$pos;
}
}
}
if($strand eq "-"){
my $pos=$start;
for(my $i=0;$i<@base;$i++){
if($base[$i] ne "-"){
$pos++;
my $real_pos = $sequence_length-1-($pos-1)+1;
$pos{$i}=$real_pos;
}
}
}
}
}
if($isRefFound == 0){
print E $content;
next;
# die "$ref not found!\nspecies name should be added before chr such as chr01 should be cattle.chr01\n";
}
my %data;
foreach my $line(@line){
my @a=split(/\s+/,$line);
my ($s,$chr,$start,$alignment_length,$strand,$sequence_length,$sequence)=@a;
$chr=~/^([^\.]+)\.(.*)/;
my $species=$1;
$chr=$2;
my @base=split(//,$sequence);
$sequence =~ tr/ATCGRYMK/TAGCYRKM/ if ($strand eq "-");
for(my $i=0;$i<@base;$i++){
next if(!exists $pos{$i});
my $pos=$pos{$i};
$data{$pos}{$species}=$base[$i];
}
}
my @output_line;
foreach my $pos(sort {$a<=>$b} keys %data){
@output_line=($ref_chr,$pos);
foreach my $species(@species){
my $base="-";
if(exists $data{$pos}{$species}){
$base = $data{$pos}{$species};
}
push @output_line,$base;
}
my $output_line=join "\t",@output_line;
print O "$output_line\n";
}
}
