IGV脚本批量截图

IGV（Integrative Genomics Viewer）是一款本地即可使用的基因组浏览器。我们常常用它来查看一些位点或者区域的变化。在一些文章的附录中也常常见到一些截取自IGV的图片。

IGV的介绍网上已经比较多了。这里就不再多说了。

一般我截图的时候是在IGV上手动点击，然后拿个截图工具一个一个截图，刚开始觉得还好，到后来需要截图的位点多了之后出现了三个问题：一就是位点太多，定位加截图比较麻烦；二就是有时候会调整一些选项导致前后截图的样式不一样；三就是截图每次截取的区域有差别，结果导致有的截图白边大，有的小，有的裁剪过头了，大小不一。后来再用ppt或者ps调整又多了一道工序就麻烦了。

手动截图的问题

在很多软件中其实都有这种批处理的命令。在IGV中有一个选项就是Run Batch Script...，就是可以执行脚本的一个命令。下面是一个脚本示例

1. 示例脚本

new
genome hg18
load  myfile.bam
snapshotDirectory mySnapshotDirectory
goto chr1:65,289,335-65,309,335
sort position
collapse
snapshot
goto chr1:113,144,120-113,164,120
sort base
collapse
snapshot sample.png
goto chr4:68,457,006-68,467,006
sort strand
collapse
snapshot
exit

上面的代码是来自IGV的官网。整个流程其实就是这样

1. 新建场景：新建 一个场景
2. 参考路径：设置参考基因组的文件路径
3. 加载文件：加载比对后得到的bam文件
4. 截图路径：设置截图的保存路径
5. 定位区域：定位到哪条染色体的哪个区域或者位置
6. 选项调整：调整一些选项（例如read排序方式、read显示方式等等）
7. 截图保存：截图并且保存（可以指定文件名，也可以不指定）
8. 退出IGV

这个过程与我们手动去点👉的步骤其实是一样的。只不过将这种点击的动作变成这种指令写到脚本中让IGV自己去读取然后执行。这样一来也就实现了自动化啦😉。下面看一下具体的步骤是怎么做的。

2. 具体使用的步骤

2.0. 查看比对情况

用IGV打开文件查看比对查看，大致了解一下情况，文件是否有效，是否比对成功等等之类的。

2.1. 写IGV的脚本

在写之前需要知道三个文件路径，

• 参考文件的路径
• 比对的bam文件的路径
• 截图保存的文件夹的路径

⏰：如果参考文件事先用IGV的load genome from file...手动加载了之后，在代码里面就可以直接输入这个基因的名字，不需要输入完整的路径了。

下面是代码以及注释，注释可以使用#，也可以使用//开头。

# 新建一个场景
new

# 设置路径（使用绝对路径，参考基因组如果之前手动加载过就可以只指定文件名）
  # 设置参考基因组的路径
genome /Volumes/HDD/Oryza.fa
  # 加载比对文件的路径（可以加载多个）
load /Volumes/HDD/Oryza.bam
  # 设置截图的文件夹路径
snapshotDirectory /Volumes/HDD/screenshot

# 定位区域
# 可以是定位到点，也可以是范围
# goto chr1:13000
goto chr1:12000-14000

# 调整选项
  # 这里将界面显示调整为read区域全显示
squish
  # 按照位置排序
sort position

# 截图
# 指定文件名的写法   snapshot 123.png
# 不制定文件名的写法 snapshot
snapshot

# 退出
exit

👉说明：在调整选项里面有两种设置比较常用。还有更多的选项，可以在文末表格中查看。

• 第一种是read显示效果，扩展、折叠、压扁。如果之前用过IGV，就会发现这个与IGV右键菜单是相似的。
  对应关系是

# 扩展
expand
# 折叠
collapse
# 压扁
squish

显示效果

• 第二种是read的排序方式，一种是按照位置排序、一种是按照碱基排序、一种是按照链的方向排序
  对应关系是

# 按照位置排序
sort position
# 按照碱基排序
sort base
# 按照链的方向排序
sort strand

上面的选项直接在IGV代码中一行写一个就可以。

2.2. 执行脚本

打开IGV>调整窗口左右的长度>然后点击Tools>Run Batch Script...>点击之前写的脚本。然后IGV会与之前我们手动点击的一样，一步一步的发生变化。最后得到截图💪。

其实为了说明，上面的我们只执行了一次截图操作，实际的应用不会只是截取一两张图片，那如果需要截取的区域较多该怎么办呢？

截图

IGV脚本比较简洁，里面没有循环之类的方法，也就是说不能写循环语句来跳转到对应的位置。每次一行一句话，按顺序执行。既然不能写循环那转换一下思路，用其他语言写的循环语句输出来生成IGV执行脚本不就行了！nice😆

3. 额外的 —— 生成IGV的执行脚本

有的时候可能需要截取的位点比较多，每次都写一句goto ...就比较麻烦，那么如何
这个时候就可以搭配其他编程语言啦。比如Perl、shell、python、java、C、Ruby等等。

3.1. 方法一：生成IGV执行的脚本

使用其他语言来生成用于IGV执行的脚本

假如我有这些位点需要截图，这些位点信息是这样的，文件名为123.txt

chr1 12000
chr1 20000
chr1 40000
chr1 70000
chr2 12000
chr2 20000
chr2 40000
chr2 70000

• 使用perl生成用来转换IGV代码的脚本Transform_IGV_batch.pl

use strict;
use warnings;
# use Getopt::Long;

# 这里为了方便演示就不使用Getopt来获得参数了。直接按照顺序来
my $reference = shift(@ARGV);  # 参考文件的路径
my $snapshot  = shift(@ARGV);  # 截图文件保存的路径
my $regionfile= shift(@ARGV);  # 需要截图的区域的文件
my @align     = @ARGV;         # 比对文件的路径(可以指定多个)

# 读取需要截图的区域的文件
my %region = ();
open my $read_region,"<",$regionfile or die $!;
while(<$read_region>){
    chomp;
    unless(m/^$/){
        my @F = split /\s+/,$_;
        push @{$region{$F[0]}},$F[1];
    }
}

print IGV_new();
print path($reference,$snapshot);
for my $load_file (@align){
    print load($load_file);
}
print preset();

for my $chr (keys %region ){
    for my $region (@{$region{$chr}}) {
        print snapshot($chr,$region);
    }
}

print IGV_exit();


sub path {
    my ($reference,$snapshot) = @_;
    my $str = "genome $reference\n";
    $str   .= "snapshotDirectory $snapshot\n";
    return $str;
}

sub load {
    my $path = shift;
    return "load $path\n";
}

sub preset {
    my @arguments = @_;
    my $str;
    if ( @arguments ){
        for my $i (@arguments){
            $str .= "$i\n";
        }
    }else{
        $str = "sort position\ncollapse\n";
    }
    return $str;
}

sub snapshot {
    my $chr    = shift;
    my $region = shift;
    my $file_name = shift || "";
    my $str = "goto ${chr}:${region}\n";
    $str   .= "snapshot $file_name\n";
    return $str;
}

sub IGV_new {
    return "new\n";
}

sub IGV_exit {
    return "exit\n";
}

使用这个脚本将上面需要截取的文件转换为IGV执行的脚本IGV_batch_script.txt

perl ./Transform_IGV_bash.pl /Volumes/HDD/Oryza.fa /Volumes/HDD/screenshot 123.txt /Volumes/HDD/Oryza.bam > IGV_batch_script.txt

然后打开IGV执行这个脚本就可以了。这个perl脚本不完善，但是基本上可以使用它来生成一个IGV的执行脚本了。里面的预设可以自己改一下。然后可以加一个getopt参数比较明确。

• 使用shell生成用来IGV执行的脚本Transform_IGV_batch.sh

#!usr/bin/bash

# 参考序列的路径
reference=$1
# 截图的路径
snapshot=$2
# 需要截图的区域的路径
regionfile=$3
# 比对的文件(可以指定多个)
n=-1
for i in "$@";
do
  ((n++))
  align[$n]=$i
done

echo "new"
echo "genome $reference"
echo "snapshotDirectory $snapshot"
for i in "${align[@]}";
do
  echo "load $i"
done

echo "sort position"
echo "collapse"

cat ${regionfile} | grep -v "^\n$" | while read IFS=$'\s' row;
do
  # regionfile是这样的格式，也可以是别的格式，那样解析的代码需要改写
  # chr1 12000
  # chr2 20000
  # chr3 40000
  chr=${row[0]}
  region=${row[1]}
  echo "goto ${chr}:${region}"
  echo "snapshot"
done

echo "exit"

使用这个脚本将上面需要截取的文件转换为IGV执行的脚本IGV_batch_script.txt

bash ./Transform_IGV_bash.sh /Volumes/HDD/Oryza.fa /Volumes/HDD/screenshot 123.txt /Volumes/HDD/Oryza.bam > IGV_batch_script.txt

下面的python脚本先占个位，以后学了再补起来。

• 使用python生成用来IGV执行的脚本Transform_IGV_batch.py

3.2. 方法二：使用接口连接到IGV

具体可以参考IGV - python这个网站给出的信息。后面我会补充起来。

4. IGV脚本支持的方法

下面是官网上给出的脚本可用方法的名字，我没有来得及对它进行翻译。

参考

• IGV - python
• Question: Automated Plotting Of Seq. Reads In A Bam File
• IGV官网