IGV腳本批量截圖

IGV（Integrative Genomics Viewer）是一款本地即可使用的基因組瀏覽器筝尾。我們常常用它來查看一些位點或者區(qū)域的變化碎赢。在一些文章的附錄中也常常見到一些截取自IGV的圖片籽孙。

IGV的介紹網(wǎng)上已經(jīng)比較多了。這里就不再多說了闷板。

一般我截圖的時候是在IGV上手動點擊逗柴，然后拿個截圖工具一個一個截圖蛹头，剛開始覺得還好，到后來需要截圖的位點多了之后出現(xiàn)了三個問題：一就是位點太多戏溺，定位加截圖比較麻煩渣蜗；二就是有時候會調(diào)整一些選項導(dǎo)致前后截圖的樣式不一樣；三就是截圖每次截取的區(qū)域有差別旷祸，結(jié)果導(dǎo)致有的截圖白邊大耕拷，有的小，有的裁剪過頭了托享，大小不一骚烧。后來再用ppt或者ps調(diào)整又多了一道工序就麻煩了。

手動截圖的問題

在很多軟件中其實都有這種批處理的命令闰围。在IGV中有一個選項就是Run Batch Script...赃绊，就是可以執(zhí)行腳本的一個命令。下面是一個腳本示例

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的官網(wǎng)羡榴。整個流程其實就是這樣

1. 新建場景：新建 一個場景
2. 參考路徑：設(shè)置參考基因組的文件路徑
3. 加載文件：加載比對后得到的bam文件
4. 截圖路徑：設(shè)置截圖的保存路徑
5. 定位區(qū)域：定位到哪條染色體的哪個區(qū)域或者位置
6. 選項調(diào)整：調(diào)整一些選項（例如read排序方式碧查、read顯示方式等等）
7. 截圖保存：截圖并且保存（可以指定文件名，也可以不指定）
8. 退出IGV

這個過程與我們手動去點??的步驟其實是一樣的炕矮。只不過將這種點擊的動作變成這種指令寫到腳本中讓IGV自己去讀取然后執(zhí)行么夫。這樣一來也就實現(xiàn)了自動化啦??。下面看一下具體的步驟是怎么做的肤视。

2. 具體使用的步驟

2.0. 查看比對情況

用IGV打開文件查看比對查看档痪，大致了解一下情況，文件是否有效邢滑，是否比對成功等等之類的腐螟。

2.1. 寫IGV的腳本

在寫之前需要知道三個文件路徑，

• 參考文件的路徑
• 比對的bam文件的路徑
• 截圖保存的文件夾的路徑

?：如果參考文件事先用IGV的load genome from file...手動加載了之后困后，在代碼里面就可以直接輸入這個基因的名字乐纸，不需要輸入完整的路徑了。

下面是代碼以及注釋摇予，注釋可以使用#汽绢，也可以使用//開頭。

# 新建一個場景
new

# 設(shè)置路徑（使用絕對路徑侧戴，參考基因組如果之前手動加載過就可以只指定文件名）
  # 設(shè)置參考基因組的路徑
genome /Volumes/HDD/Oryza.fa
  # 加載比對文件的路徑（可以加載多個）
load /Volumes/HDD/Oryza.bam
  # 設(shè)置截圖的文件夾路徑
snapshotDirectory /Volumes/HDD/screenshot

# 定位區(qū)域
# 可以是定位到點宁昭，也可以是范圍
# goto chr1:13000
goto chr1:12000-14000

# 調(diào)整選項
  # 這里將界面顯示調(diào)整為read區(qū)域全顯示
squish
  # 按照位置排序
sort position

# 截圖
# 指定文件名的寫法   snapshot 123.png
# 不制定文件名的寫法 snapshot
snapshot

# 退出
exit

??說明：在調(diào)整選項里面有兩種設(shè)置比較常用跌宛。還有更多的選項，可以在文末表格中查看积仗。

• 第一種是read顯示效果疆拘，擴展、折疊寂曹、壓扁哎迄。如果之前用過IGV，就會發(fā)現(xiàn)這個與IGV右鍵菜單是相似的隆圆。
  對應(yīng)關(guān)系是

# 擴展
expand
# 折疊
collapse
# 壓扁
squish

顯示效果

• 第二種是read的排序方式漱挚，一種是按照位置排序、一種是按照堿基排序匾灶、一種是按照鏈的方向排序
  對應(yīng)關(guān)系是

# 按照位置排序
sort position
# 按照堿基排序
sort base
# 按照鏈的方向排序
sort strand

上面的選項直接在IGV代碼中一行寫一個就可以棱烂。

2.2. 執(zhí)行腳本

打開IGV>調(diào)整窗口左右的長度>然后點擊Tools>Run Batch Script...>點擊之前寫的腳本租漂。然后IGV會與之前我們手動點擊的一樣阶女，一步一步的發(fā)生變化。最后得到截圖??哩治。

其實為了說明秃踩，上面的我們只執(zhí)行了一次截圖操作，實際的應(yīng)用不會只是截取一兩張圖片业筏，那如果需要截取的區(qū)域較多該怎么辦呢憔杨？

截圖

IGV腳本比較簡潔，里面沒有循環(huán)之類的方法蒜胖，也就是說不能寫循環(huán)語句來跳轉(zhuǎn)到對應(yīng)的位置消别。每次一行一句話，按順序執(zhí)行台谢。既然不能寫循環(huán)那轉(zhuǎn)換一下思路寻狂，用其他語言寫的循環(huán)語句輸出來生成IGV執(zhí)行腳本不就行了！nice??

3. 額外的 —— 生成IGV的執(zhí)行腳本

有的時候可能需要截取的位點比較多朋沮，每次都寫一句goto ...就比較麻煩蛇券，那么如何
這個時候就可以搭配其他編程語言啦。比如Perl樊拓、shell纠亚、python、java筋夏、C蒂胞、Ruby等等。

3.1. 方法一：生成IGV執(zhí)行的腳本

使用其他語言來生成用于IGV執(zhí)行的腳本

假如我有這些位點需要截圖条篷，這些位點信息是這樣的骗随，文件名為123.txt

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

• 使用perl生成用來轉(zhuǎn)換IGV代碼的腳本Transform_IGV_batch.pl

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

# 這里為了方便演示就不使用Getopt來獲得參數(shù)了岳瞭。直接按照順序來
my $reference = shift(@ARGV);  # 參考文件的路徑
my $snapshot  = shift(@ARGV);  # 截圖文件保存的路徑
my $regionfile= shift(@ARGV);  # 需要截圖的區(qū)域的文件
my @align     = @ARGV;         # 比對文件的路徑(可以指定多個)

# 讀取需要截圖的區(qū)域的文件
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";
}

使用這個腳本將上面需要截取的文件轉(zhuǎn)換為IGV執(zhí)行的腳本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執(zhí)行這個腳本就可以了。這個perl腳本不完善蚊锹，但是基本上可以使用它來生成一個IGV的執(zhí)行腳本了瞳筏。里面的預(yù)設(shè)可以自己改一下。然后可以加一個getopt參數(shù)比較明確牡昆。

• 使用shell生成用來IGV執(zhí)行的腳本Transform_IGV_batch.sh

#!usr/bin/bash

# 參考序列的路徑
reference=$1
# 截圖的路徑
snapshot=$2
# 需要截圖的區(qū)域的路徑
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"

使用這個腳本將上面需要截取的文件轉(zhuǎn)換為IGV執(zhí)行的腳本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腳本先占個位丢烘，以后學(xué)了再補起來柱宦。

• 使用python生成用來IGV執(zhí)行的腳本Transform_IGV_batch.py

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

具體可以參考IGV - python這個網(wǎng)站給出的信息。后面我會補充起來播瞳。

4. IGV腳本支持的方法

下面是官網(wǎng)上給出的腳本可用方法的名字掸刊，我沒有來得及對它進行翻譯。

參考

• IGV - python
• Question: Automated Plotting Of Seq. Reads In A Bam File
• IGV官網(wǎng)