批量IGV截图【直播】我的基因组83
把我的全基因组重测续数据bam文件载入到IGV看了几个基因，发现有一些基因比对情况非常诡异，各种色块，各种缺口，让我不忍直视，搞得像是个破损的基因组，也查了查那些基因，主要是一些家族性基因，太长的基因，或者复杂度非常低的基因。所以我就想看看其它基因如何，本来是准备一个个查询，截屏的，但是所有批量操作都应该是可以编程解决的，就搜索了一下，的确找到了解决方案。成功的截屏了50000张IGV图片。
首先从gencode数据库里面可以下载所有的gtf文件
下载代码是：
mkdir -p ~/reference/gtf/gencodecd ~/reference/gtf/gencode## https://www.gencodegenes.org/releases/current.htmlwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.2wayconspseudos.gtf.gzwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.long_noncoding_RNAs.gtf.gzwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.polyAs.gtf.gzwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gtf.gz## https://www.gencodegenes.org/releases/25lift37.htmlwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.annotation.gtf.gzwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.HGNC.gzwget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.EntrezGene.gzwget 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.positionzcat 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.positionzcat 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.positionzcat 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.positionzcat 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.positionzcat 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.positionchr1 69091 70008 OR4F5chr1 367640 368634 OR4F29chr1 621096 622034 OR4F16chr1 859308 879961 SAMD11chr1 879584 894689 NOC2Lchr1 895967 901095 KLHL17chr1 901877 911245 PLEKHN1chr1 910584 917473 PERM1chr1 934342 935552 HES4chr1 936518 949921 ISG15
批量IGV截图脚本
这里我想把我的全基因组重测续数据载入到IGV并且根据所有的基因查看并且截图保存。其中IGV脚本很简单，格式如下：
goto chr1:25158588-25161609snapshot chr1_25158688_25161509_slop100.pnggoto chr1:26489713-26490238snapshot 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.txtgoto chr1:7745384-7775384snapshot CAMTA1_31.pnggoto chr1:7775384-7805384snapshot CAMTA1_32.pnggoto chr1:7805384-7835384snapshot CAMTA1_33.pnggoto chr1:7831329-7841492snapshot VAMP3.pnggoto chr1:7844380-7874380snapshot PER3_1.pnggoto chr1:7874380-7904380snapshot PER3_2.pnggoto chr1:7904380-7934380snapshot PER3_3.pnggoto chr1:7903143-7913572snapshot UTS2.pnggoto chr1:7975954-8003225snapshot TNFRSF9.pnggoto chr1:8014351-8044351snapshot PARK7_1.pnggoto chr1:8044351-8074351snapshot PARK7_2.pnggoto chr1:8064464-8086368snapshot ERRFI1.png
比如CAMTA1这个基因长约1M，所以会被切割成33个片段才出图。
这个IGV脚本运行指导流程是：
运行结果如下：