批量IGV截图【直播】我的基因组83

查看原文

其他

批量IGV截图【直播】我的基因组83

原创 2017-07-18 曾健明 生信技能树

把我的全基因组重测续数据bam文件载入到IGV看了几个基因，发现有一些基因比对情况非常诡异，各种色块，各种缺口，让我不忍直视，搞得像是个破损的基因组，也查了查那些基因，主要是一些家族性基因，太长的基因，或者复杂度非常低的基因。所以我就想看看其它基因如何，本来是准备一个个查询，截屏的，但是所有批量操作都应该是可以编程解决的，就搜索了一下，的确找到了解决方案。成功的截屏了50000张IGV图片。

首先从gencode数据库里面可以下载所有的gtf文件

下载代码是：

mkdir -p ~/reference/gtf/gencode
cd ~/reference/gtf/gencode
## https://www.gencodegenes.org/releases/current.html
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.2wayconspseudos.gtf.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.long_noncoding_RNAs.gtf.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.polyAs.gtf.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gtf.gz
## https://www.gencodegenes.org/releases/25lift37.html
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.annotation.gtf.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.HGNC.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.EntrezGene.gz
wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.RefSeq.gz

然后写脚本得到基因的染色体还有起始终止坐标

代码是：

zcat gencode.v25.long_noncoding_RNAs.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >lncRNA.hg38.position
zcat gencode.v25.2wayconspseudos.gtf.gz |perl -alne '{next unless $F[2] eq "transcript" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >pseudos.hg38.position
zcat gencode.v25.annotation.gtf.gz| grep protein_coding |perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >protein_coding.hg38.position
zcat gencode.v25.annotation.gtf.gz|perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >allGene.hg38.position
zcat gencode.v25lift37.annotation.gtf.gz | grep protein_coding |perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >protein_coding.hg19.position
zcat gencode.v25lift37.annotation.gtf.gz | perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >allGene.hg19.position

PS：这里面有一个小问题，gencode里面的数据有着HAVANA和ENSEMBL的区别，尤其是在hg38里面，需要区别对待！

查看基因的坐标信息bed文件如

head ~/reference/gtf/gencode/protein_coding.hg19.position
chr1 69091 70008 OR4F5
chr1 367640 368634 OR4F29
chr1 621096 622034 OR4F16
chr1 859308 879961 SAMD11
chr1 879584 894689 NOC2L
chr1 895967 901095 KLHL17
chr1 901877 911245 PLEKHN1
chr1 910584 917473 PERM1
chr1 934342 935552 HES4
chr1 936518 949921 ISG15

批量IGV截图脚本

这里我想把我的全基因组重测续数据载入到IGV并且根据所有的基因查看并且截图保存。其中IGV脚本很简单，格式如下：

goto chr1:25158588-25161609
snapshot chr1_25158688_25161509_slop100.png
goto chr1:26489713-26490238
snapshot chr1_26489813_26490138_slop100.png

参考：https://software.broadinstitute.org/software/igv/batch

就是进入基因的区域，然后截图保存即可。但是要注意IGV对bam文件的分辨率上线是37kb。但是人类的近2万个编码蛋白的基因里面有接近一半都超过了这个上线，所以这些基因都需要截多张图！

制作这个IGV脚本的代码如下：

cat ~/reference/gtf/gencode/protein_coding.hg19.position| \
perl -alne '{if (($F[2]-$F[1])<30000){print "goto $F[0]:$F[1]-$F[2]\nsnapshot $F[3].png"} \
else{$n=int(($F[2]-$F[1])/30000)+1;foreach (1..$n){$start=$F[1]+($_-1)*30000;$end=$start+30000;print "goto $F[0]:$start-$end\nsnapshot $F[3]_$_.png" }}}' \
>igv_all_gene_snapshot.txt
goto chr1:7745384-7775384
snapshot CAMTA1_31.png
goto chr1:7775384-7805384
snapshot CAMTA1_32.png
goto chr1:7805384-7835384
snapshot CAMTA1_33.png
goto chr1:7831329-7841492
snapshot VAMP3.png
goto chr1:7844380-7874380
snapshot PER3_1.png
goto chr1:7874380-7904380
snapshot PER3_2.png
goto chr1:7904380-7934380
snapshot PER3_3.png
goto chr1:7903143-7913572
snapshot UTS2.png
goto chr1:7975954-8003225
snapshot TNFRSF9.png
goto chr1:8014351-8044351
snapshot PARK7_1.png
goto chr1:8044351-8074351
snapshot PARK7_2.png
goto chr1:8064464-8086368
snapshot ERRFI1.png

比如CAMTA1这个基因长约1M，所以会被切割成33个片段才出图。

这个IGV脚本运行指导流程是：

运行结果如下：

直播我的基因组分析-目录-持续更新

【直播】我的基因组81：看看我的vcf文件的vaf分布情况

【直播】我的基因组82：如何对maf格式的突变文件统计vaf

李尚福、魏凤和双双被拿下，与美国一份报告是否有关？

这样的洞庭湖决堤，实在让人同情不起来

有的人走了，却永远活着

圈内疯传某谣言

不要放过这些人渣

批量IGV截图【直播】我的基因组83

首先从gencode数据库里面可以下载所有的gtf文件

然后写脚本得到基因的染色体还有起始终止坐标

查看基因的坐标信息bed文件如

批量IGV截图脚本

您可能也对以下帖子感兴趣

李尚福、魏凤和双双被拿下，与美国一份报告是否有关？

这样的洞庭湖决堤，实在让人同情不起来

有的人走了，却永远活着

圈内疯传某谣言

不要放过这些人渣

生成图片，分享到微信朋友圈

批量IGV截图【直播】我的基因组83

首先从gencode数据库里面可以下载所有的gtf文件

然后写脚本得到基因的染色体还有起始终止坐标

查看基因的坐标信息bed文件如

批量IGV截图脚本

您可能也对以下帖子感兴趣