IGV脚本批量截图
2019-01-02 本文已影响11人
白菜代码小推车
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的官网。整个流程其实就是这样
- 新建场景：新建 一个场景
- 参考路径：设置参考基因组的文件路径
- 加载文件：加载比对后得到的bam文件
- 截图路径：设置截图的保存路径
- 定位区域：定位到哪条染色体的哪个区域或者位置
- 选项调整：调整一些选项（例如read排序方式、read显示方式等等）
- 截图保存：截图并且保存（可以指定文件名，也可以不指定）
- 退出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脚本支持的方法
下面是官网上给出的脚本可用方法的名字，我没有来得及对它进行翻译。
| Command | Description |
|---|---|
| new | Create a new session. Unloads all tracks except the default genome annotations. |
| load file | Loads data or session files. Specify a comma-delimited list of full paths or URLs.Note: for Google Genomics readgroup sets the id is the "file", specify format=ga4gh (version 2.3.80 and greater only). For exampleload CMvnhpKTFhCjz9_25e_lCw format=ga4gh To explicitly specify a path to an index file use the optional "index=" parameter. load foo.bam index=bar.bai |
| collapse trackName | Collapses a given trackName. trackName is optional, however, and if it is not supplied all tracks are collapsed. |
| echo | Writes "echo" back to the response. (Primarily for testing) |
| exit | Exit (close) the IGV application. |
| expand trackName | Expands a given trackName. trackName is optional, however, and if it is not supplied all tracks are expanded. |
| genome genomeIdOrPath | Selects a genome by id, or loads a genome (or indexed fasta) from the supplied path. |
| goto locus or listOfLoci | Scrolls to a single locus or a space-delimited list of loci. If a list is provided, these loci will be displayed in a split screen view. Use any syntax that is valid in the IGV search box. |
| goto all | Scrolls to a whole genome view. |
| region chr start end | Defines a region of interest bounded by the two loci (e.g., region chr1 100 200). |
| maxPanelHeight height | Sets the number of vertical pixels (height) of each panel to include in image. Images created from a port command or batch script are not limited to the data visible on the screen. Stated another way, images can include the entire panel not just the portion visible in the scrollable screen area. The default value for this setting is 1000, increase it to see more data, decrease it to create smaller images. |
| setLogScale(true | false) | |
| setSleepInterval ms | Sets a delay (sleep) time in milliseconds. The sleep interval is invoked between successive commands. |
| snapshotDirectory path | Sets the directory in which to write images. |
| snapshot filename | Saves a snapshot of the IGV window to an image file. If filename is omitted, writes a PNG file with a filename generated based on the locus. If filename is specified, the filenameextension determines the image file format, which must be .png, .jpg, or .svg. |
| sort option locus | Sorts an alignment track by the specified option. Recognized values for the option parameter are: base, position, strand, quality, sample, readGroup, AMPLIFICATION, DELETION, EXPRESSION, SCORE, and MUTATION_COUNT. The locus option can define a single position, or a range. If absent sorting will be perfomed based on the region in view, or the center position of the region in view, depending on the option. |
| squish trackName | Squish a given trackName. trackName is optional, and if it is not supplied all annotation tracks are squished. |
| viewaspairs trackName | Set the display mode for an alignment track to "View as pairs". trackName is optional. |
| preference key value | Temporarily set the preference named key to the specified value. This preference only lasts until IGV is shut down. |
