搜索
查看: 2857|回复: 0

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

[复制链接]

23

主题

37

帖子

370

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
370
发表于 2017-7-19 16:09:05 | 显示全部楼层 |阅读模式
批量IGV截图【直播】我的基因组83  原文链接

把我的全基因组重测续数据bam文件载入到IGV看了几个基因,发现有一些基因比对情况非常诡异,各种色块,各种缺口,让我不忍直视,搞得像是个破损的基因组,也查了查那些基因,主要是一些家族性基因,太长的基因,或者复杂度非常低的基因。所以我就想看看其它基因如何,本来是准备一个个查询,截屏的,但是所有批量操作都应该是可以编程解决的,就搜索了一下,的确找到了解决方案。成功的截屏了50000张IGV图片。首先从gencode数据库里面可以下载所有的gtf文件下载代码:
[AppleScript] 纯文本查看 复制代码
mkdir -p ~/reference/gtf/gencode
cd  ~/reference/gtf/gencode
## [url=https://www.gencodegenes.org/releases/current.html]https://www.gencodegenes.org/releases/current.html[/url]
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.2wayconspseudos.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... yconspseudos.gtf.gz[/url]
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.long_noncoding_RNAs.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... ncoding_RNAs.gtf.gz[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.polyAs.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... e.v25.polyAs.gtf.gz[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... 5.annotation.gtf.gz[/url] 
## [url=https://www.gencodegenes.org/releases/25lift37.html]https://www.gencodegenes.org/releases/25lift37.html[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.annotation.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... 7.annotation.gtf.gz[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.HGNC.gz]ftp://ftp.sanger.ac.uk/pub/genco ... 37.metadata.HGNC.gz[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.EntrezGene.gz]ftp://ftp.sanger.ac.uk/pub/genco ... adata.EntrezGene.gz[/url] 
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.RefSeq.gz]ftp://ftp.sanger.ac.uk/pub/genco ... .metadata.RefSeq.gz[/url] 
然后写脚本得到基因的染色体还有起始终止坐标
[AppleScript] 纯文本查看 复制代码
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文件如
[AppleScript] 纯文本查看 复制代码
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脚本很简单,格式如下:
[AppleScript] 纯文本查看 复制代码
goto chr1:25158588-25161609
snapshot chr1_25158688_25161509_slop100.png
goto chr1:26489713-26490238
snapshot chr1_26489813_26490138_slop100.png
就是进入基因的区域,然后截图保存即可。但是要注意IGV对bam文件的分辨率上线是37kb。但是人类的近2万个编码蛋白的基因里面有接近一半都超过了这个上线,所以这些基因都需要截多张图!
制作这个IGV脚本的代码如下:
[AppleScript] 纯文本查看 复制代码
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脚本运行指导流程是:

运行结果如下:










上一篇:biplot 画主成分分析聚类图小例子
下一篇:亚马逊云教程1:云是什么,创建并连接EC2
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:34 , Processed in 0.037598 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.