Calculate phastCon Score for a gene —- 計算基因的phastCon平均分，判斷基因保守型

PhastCon socre is the score from 0 to 1 to show the conservation level.

A score showing the posterior probability that phastCons’s phylogenetic hidden Markov model (HMM) is in its most conserved state at that base position.

The phastCons scores represent probabilities of negative selection and range between 0 and 1.

Short highly-conserved regions and long moderately conserved regions can both obtain high scores.

也就是說如果某個位點或者一段序列的phastCon分值高的話，表示保守型較高。

1，下載phast score的wigFix格式文件(以19號染色體為測試)

wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phastCons46way/placentalMammals/chr19.phastCons46way.placental.wigFix.gz

gunzip -k chr19.phastCons46way.placental.wigFix.gz

head chr19.phastCons46way.placental.wigFix

fixedStep chrom=chr19 start=60001 step=1

0.412

0.407

0.394

0.393

0.381

0.381

0.369

0.345

0.291

wigFix格式是 fixed 的 Wiggle 格式，根據上面的前幾行，應該能夠看出該文件表示19號染色體，從60001位置開始，每一步step(base)的phastCon分值。

利用bedops，將wig轉成bed文件

wig2bed chr19.phastCons46way.placental.wigFix > chr19.phastCons46way.placental.bed

如果需要合并的話，運行下面的命令即可，這里以chr19染色體為例，就不合并了。

bedops –everything chr*.bed > vertebrate.phastCons46.bed

將gtf文件分割成已基因為單位的gtf文件

gtf文件中只有chr19染色體的基因注釋信息。

mkdir test

awk -F 't' '{split($9,a,";"); split(a[1],a," "); gsub(""", "", a[2]); x="test/"a[2]".gtf"; print >> ("test/"a[2]".gtf"); close(x);}' Homo_sapiens.GRCh38.85.gtf && rm test/.gtf

將gft文件轉成bed文件

rm test/*bed

for fn in `ls test/*gtf`; do

gene=`echo $fn | cut -d . -f 1`

gtf2bed $fn | sort -k1,1 -k2,2n >> $gene.bed

rm $fn

done

gtf2bed是前面文章中提到的，只要能將gtf轉成bed就可。

利用bedtools工具計算每個序列的平均phastCon score

# 最新版的bedtools有split選項

for fn in `ls test/*.bed`;do

bedtools map -split -sorted -a $fn -b chr19.phastCons46way.placental.bed -o mean | awk -F 't' '{printf ("%st%4fn", $1,$NF)}'

done

bedtools官網 bedtools.readthedocs.io ，mean表示計算平均分。結果如下，平均分越高越保守。

chr19 531711 542092 ENST00000215574 0.388707

chr19 531768 536715 ENST00000586788 0.317417

chr19 532031 542092 ENST00000586283 0.431241

chr19 532099 541585 ENST00000607527 0.643244

chr19 535836 542087 ENST00000606065 0.404337

chr19 536952 542042 ENST00000606400 0.288264

chr19 541216 541661 ENST00000593036 0.434159

chr19 489175 505342 ENST00000587541 0.032945

chr19 496453 505207 ENST00000613880 0.186624

chr19 496453 505207 ENST00000617201 0.194969

參考：

https://www.biostars.org/p/150152/

http://ccg.vital-it.ch/mga/mm9/phastcons/phastcons.html

#####################################################################

#版權所有 轉載請告知 版權歸作者所有 如有侵權 一經發現 必將追究其法律責任

#Author: Jason

####################################################################

相關

總結

以上是生活随笔為你收集整理的html计算平均分,Calculate phastCon Score for a gene —- 计算基因的phastCon平均分，判断基因保守型...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。