Calculate phastCon Score for a gene —- 計(jì)算基因的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.

也就是說(shuō)如果某個(gè)位點(diǎn)或者一段序列的phastCon分值高的話，表示保守型較高。

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

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 格式，根據(jù)上面的前幾行，應(yīng)該能夠看出該文件表示19號(hào)染色體，從60001位置開(kāi)始，每一步step(base)的phastCon分值。

利用bedops，將wig轉(zhuǎn)成bed文件

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

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

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

將gtf文件分割成已基因?yàn)閱挝坏膅tf文件

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文件轉(zhuǎn)成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轉(zhuǎn)成bed就可。

利用bedtools工具計(jì)算每個(gè)序列的平均phastCon score

# 最新版的bedtools有split選項(xiàng)

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官網(wǎng) bedtools.readthedocs.io ，mean表示計(jì)算平均分。結(jié)果如下，平均分越高越保守。

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

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

#版權(quán)所有 轉(zhuǎn)載請(qǐng)告知 版權(quán)歸作者所有 如有侵權(quán) 一經(jīng)發(fā)現(xiàn) 必將追究其法律責(zé)任

#Author: Jason

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

相關(guān)

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。