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

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

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

下载代码是：

  1. mkdir -p ~/reference/gtf/gencode

  2. cd  ~/reference/gtf/gencode

  3. ## https://www.gencodegenes.org/releases/current.html

  4. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.2wayconspseudos.gtf.gz

  5. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.long_noncoding_RNAs.gtf.gz

  6. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.polyAs.gtf.gz

  7. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gtf.gz

  8. ## https://www.gencodegenes.org/releases/25lift37.html

  9. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.annotation.gtf.gz

  10. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.HGNC.gz

  11. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.EntrezGene.gz

  12. wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/GRCh37_mapping/gencode.v25lift37.metadata.RefSeq.gz

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

代码是：

  1. 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

  2. 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

  3. 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

  4. 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

  5. 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

  6. 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文件如

  1. head ~/reference/gtf/gencode/protein_coding.hg19.position

  2. chr1    69091   70008   OR4F5

  3. chr1    367640  368634  OR4F29

  4. chr1    621096  622034  OR4F16

  5. chr1    859308  879961  SAMD11

  6. chr1    879584  894689  NOC2L

  7. chr1    895967  901095  KLHL17

  8. chr1    901877  911245  PLEKHN1

  9. chr1    910584  917473  PERM1

  10. chr1    934342  935552  HES4

  11. chr1    936518  949921  ISG15

批量IGV截图脚本

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

  1. goto chr1:25158588-25161609

  2. snapshot chr1_25158688_25161509_slop100.png

  3. goto chr1:26489713-26490238

  4. snapshot chr1_26489813_26490138_slop100.png

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

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

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

  1. cat  ~/reference/gtf/gencode/protein_coding.hg19.position| \

  2. perl -alne '{if (($F[2]-$F[1])<30000){print "goto $F[0]:$F[1]-$F[2]\nsnapshot $F[3].png"} \

  3. 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" }}}' \

  4. >igv_all_gene_snapshot.txt

  5. goto    chr1:7745384-7775384

  6. snapshot        CAMTA1_31.png

  7. goto    chr1:7775384-7805384

  8. snapshot        CAMTA1_32.png

  9. goto    chr1:7805384-7835384

  10. snapshot        CAMTA1_33.png

  11. goto    chr1:7831329-7841492

  12. snapshot        VAMP3.png

  13. goto    chr1:7844380-7874380

  14. snapshot        PER3_1.png

  15. goto    chr1:7874380-7904380

  16. snapshot        PER3_2.png

  17. goto    chr1:7904380-7934380

  18. snapshot        PER3_3.png

  19. goto    chr1:7903143-7913572

  20. snapshot        UTS2.png

  21. goto    chr1:7975954-8003225

  22. snapshot        TNFRSF9.png

  23. goto    chr1:8014351-8044351

  24. snapshot        PARK7_1.png

  25. goto    chr1:8044351-8074351

  26. snapshot        PARK7_2.png

  27. goto    chr1:8064464-8086368

  28. snapshot        ERRFI1.png

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

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

运行结果如下：

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

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

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

赞 (0)

相关推荐