实用小脚本

发表于 | 更新于 | 分类于

Perl

1. 拆分FASTA文件

将FASTA文件按照用户指定的序列条数进行拆分，即每n条序列写到一个文件中，或者按照指定的输出文件数进行拆分，则每个文件中包含的序列数相同

脚本名：splitFasta.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;
use POSIX;

# 帮助文档
=head1 Description

	This script is used to split fasta file, which is too large with thosands of sequence

=head1 Usage

	$0 -i <input> -o <output_dir> [-n <seq_num_per_file>] [-m <output_file_num>]
	
=head1 Parameters

	-i	[str]	Input raw fasta file
	-o	[str]	Output file to which directory
	-n	[int]	Sequence number per file, alternate chose paramerter "-n" or "-m", if set "-n" and "-m" at the same time, only take "-n" parameter
	-m	[int]	Output file number (default:100)
=cut

my ($input,$output_dir,$seq_num,$file_num);
GetOptions(
	"i:s"=>\$input,
	"o:s"=>\$output_dir,
	"n:i"=>\$seq_num,
	"m:i"=>\$file_num
	);

die `pod2text $0` if ((!$input) or (!$output_dir));

# 设置每个文件的序列条数
if(!defined($seq_num)){
	if(!defined($file_num)){
		$file_num=100;
		my $total_seq_num=`awk 'BEGIN{n=0} /^>/{n++} END{print n}' $input`;
		chomp $total_seq_num;
		$seq_num=ceil($total_seq_num/$file_num);
	}else{
		my $total_seq_num=`awk 'BEGIN{n=0} /^>/{n++} END{print n}' $input`;
		chomp $total_seq_num;
		$seq_num=ceil($total_seq_num/$file_num);
	}
}

open IN,"<$input" or die "Cann't open $input\n";

my $n_seq=0;	# 该变量用于记录当前扫描到的序列数
my $n_file=1;	# 该变量用于记录当前真正写入的文件的计数
my $input_base=`basename $input`;
chomp $input_base;

open OUT,">$output_dir/${input_base}_${n_file}" or die "Cann't create $output_dir/${input_base}_${n_file}\n";

while(<IN>){
	next if (/^\s+$/);	# 跳过空行
	chomp;
	if (/^>/){
		$n_seq++;
		# 判断目前已经扫描到的序列数，若大于设定的split的序列数，则创建新文件
		if ($n_seq>$seq_num){
			$n_seq=1;
			$n_file++;
			close OUT;
			open OUT,">$output_dir/${input_base}_${n_file}" or die "Cann't create $output_dir/${input_base}_${n_file}\n";
			print OUT "$_\n";
		}else{
			print OUT "$_\n";
		}
	}else{
		print OUT "$_\n";
	}
}

close IN;

执行方法： 

1
$ perl splitFasta.pl -i <input> -o <output_dir> [-n <seq_num_per_file>] [-m <output_file_num>]

具体的参数使用说明可以执行 perl splitFasta.pl，查看脚本使用文档

2. 从双端FASTQ文件中抽取指定数据量(bp)的序列

按照用户指定的数据量，即多少bp，从原始的双端FASTQ文件中随机抽取序列，要求双端FASTQ文件中的序列必须配对

脚本名：extractFastq.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;

# 脚本帮助文档
=head1 Description

	Thise script is used to extract a number of fastq recorders from original input fastq file

=head1 Usage

	$0 -n <totalReads> -e <base-pairs to extract> -l <readLength> -1 <input1.fastq> -2 <input2.fastq> [-o <outdir>]

=head1 Parameters

	-n	[int]	Total reads number of input.fastq
	-e	[int]	Base-pairs per end you want to extract
	-l	[int]	The reads length
	-1	[str]	Input 1st-end fastq file
	-2	[str]	Input 2nd-end fastq file
	-o	[str]	Output directory [default: current folder]

=cut

my ($totalReads,$bpNum,$length,$Input1,$Input2,$Outdir);
GetOptions(
	"n:i"=>\$totalReads,
	"e:i"=>\$bpNum,
	"l:i"=>\$length,
	"1:s"=>\$Input1,
	"2:s"=>\$Input2,
	"o:s"=>\$Outdir
	);

$Outdir=`pwd` unless (defined($Outdir));
die `pod2text $0` if ((!$totalReads) or (!$bpNum)) or (!$length) or (!$Input1) or (!$Input2);

open FQ1,"<$Input1" or die "$!\n";
open FQ2,"<$Input2" or die "$!\n";

my $Input1_basename=`basename $Input1`;
chomp $Input1_basename;
my $Input2_basename=`basename $Input2`;
chomp $Input2_basename;
open OUT1,">$Outdir/${Input1_basename}.extract" or die "$!\n";
open OUT2,">$Outdir/${Input2_basename}.extract" or die "$!\n";

my $readsRemain=$bpNum/$length;

my $remainCount=0;

while(! eof($FQ1)){
	# 读入1st-end fastq 文件的四行
	my $fq1_1=<FQ1>;
	my $fq1_2=<FQ1>;
	my $fq1_3=<FQ1>;
	my $fq1_4=<FQ1>;
	chomp($fq1_1,$fq1_2,$fq1_3,$fq1_4);
	# 读入2nd-end fastq 文件的四行
	my $fq2_1=<FQ2>;
	my $fq2_2=<FQ2>;
	my $fq2_3=<FQ2>;
	my $fq2_4=<FQ2>;
	chomp($fq2_1,$fq2_2,$fq2_3,$fq2_4);

	# 随机抽取
	if (rand()<$readsRemain/$totalReads){
		$remainCount++;
		print OUT1 "$fq1_1\n$fq1_2\n$fq1_3\n$fq1_4\n";
		print OUT2 "$fq2_1\n$fq2_2\n$fq2_3\n$fq2_4\n";
	}
}

print "Total reads: $totalReads\n";
print "Theorical remained reads: $readsRemain\n";
print "Practical remained reads: $remainCount\n";

close FQ1;
close FQ2;
close OUT1;
close OUT2;

执行方法：

1
$ perl extractFastq.pl -n <totalReads> -e <base-pairs to extract> -l <readLength> -1 <input1.fastq> -2 <input2.fastq> [-o <outdir>]

具体的参数使用说明可以执行 perl extractFastq.pl，查看脚本使用文档

3. 双端FASTQ文件进行双端配对

脚本名：PairsMate.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;

# 脚本帮助文档
=head1 Description

Thise script is used to match the pair-end reads from one sample

=head1 Usage

$0 -1 <input1.fastq> -2 <input2.fastq> [-o <outdir>]

=head1 Parameters

-1	[str]	Input 1st-end fastq file
-2	[str]	Input 2nd-end fastq file
-o	[str]	Output directory [default: current folder]

=cut

my ($Input1,$Input2,$Outdir);
GetOptions(
"1:s"=>\$Input1,
"2:s"=>\$Input2,
"o:s"=>\$Outdir
);

$Outdir=`pwd` unless (defined($Outdir));
die `pod2text $0` if (!$Input1) or (!$Input2);

open FQ1,"<$Input1" or die "Cann't open $Input1\n";
open FQ2,"<$Input2" or die "Cann't open $Input2\n";

my $Input1_basename=`basename $Input1`;
chomp $Input1_basename;
my $Input2_basename=`basename $Input2`;
chomp $Input2_basename;
open OUT1,">$Outdir/${Input1_basename}.match" or die "Cann't create $Outdir/${Input1_basename}.match\n";
open OUT2,">$Outdir/${Input2_basename}.match" or die "Cann't create $Outdir/${Input2_basename}.match\n";

# 载入两个fastq文件，保存成哈希
my (%fq1_seq,%fq1_qua,%fq2_seq,%fq2_qua);
while (!eof(FQ1)){
        my ($id,)=split(/\s/,<FQ1>);
        $fq1_seq{$id}=<FQ1>;
        <FQ1>;
        $fq1_qua{$id}=<FQ1>;
        chomp ($fq1_seq{$id},$fq1_qua{$id});
}
while (!eof(FQ2)){
        my ($id,)=split(/\s/,<FQ2>);
        $fq2_seq{$id}=<FQ2>;
        <FQ2>;
        $fq2_qua{$id}=<FQ2>;
        chomp ($fq2_seq{$id},$fq2_qua{$id});
}
close FQ1;
close FQ2;

# 双端配对
foreach my $key (sort keys %fq1_seq){
        if(defined($fq2_seq{$key})){
                print OUT1 "$key\n$fq1_seq{$key}\n+\n$fq1_qua{$key}\n";
                print OUT2 "$key\n$fq2_seq{$key}\n+\n$fq2_qua{$key}\n";
        }
}
close OUT1;
close OUT2;

执行方法：

1
$ perl PairsMate.pl -1 <input1.fastq> -2 <input2.fastq> [-o <outdir>]

具体的参数使用说明可以执行 perl PairsMate.pl，查看脚本使用文档

4. 根据序列Id提取FASTA序列

提供FASTA文件，和包含序列Id的文件（有多列，用制表符隔开，其中一列为序列Id），提取该序列Id对应的FASTA序列

脚本名：extractSeqFromFasta.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#!/usr/bin/perl
use strict;
use warnings;

my geneList = $ARGV[0];
my fasta = $ARGV[1];
my out = $ARGV[2];

my %Hash_fasta;
my $seqId;

# 读入fasta文件，存为哈希
open FA,"<$fasta" or die "$!\n";
while(<FA>){
	chomp;
	next if(/^\s?$/);
	if(/^>(.+)$/){
		$seqId = $1;
	}else{
		$Hash_fasta{$seqId} .= $_;
	}
}
close(FA);

# 逐行读入geneList文件，并从上一步的哈希中将该序列提取出来
open LIST,"<$geneList" or die "$!\n";
open OUT,">$out" or die "$!\n";
while(<LIST>){
	chomp;
	next if(/^\s?$/);
	my $gene = $_;
	if($Hash_fasta{$gene}){
		print OUT ">$gene\n$Hash_fasta{$gene}\n";
	}
}
close(LIST);
close(OUT);

使用方法：

1
$ perl extractSeqFromFasta.pl <in.fasta> <gene list> <out.fasta>

5. 统计FastQC输出

统计FastQC结果，注意脚本需要在FastQC输出结果所在目录下运行。脚本的统计项目及输出如下：

Sample Total Reads GC Content Q20 Q30
NP006_RRS05401_1.clean.fq.extract_fastqc 39998924 49 0.999430935 0.946251454
NP006_RRS05401_2.clean.fq.extract_fastqc 39998924 49 0.986292131 0.883245809
NP007_RRS05402_1.clean.fq.extract_fastqc 40003323 49 0.999429311 0.946197936

脚本名：fastqc_stat.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
opendir (DIR, "./") or die "can't open the directory!";
@dir = readdir DIR;
foreach $file ( sort  @dir) 
{

# 跳过不需要的文件/文件夹，留下需要的文件夹
next unless -d $file;
next if $file eq '.';
next if $file eq '..';

# 提取total reads
$total_reads=  `grep '^Total' ./$file/fastqc_data.txt`;
$total_reads=(split(/\s+/,$total_reads))[2];
# 提取%GC
$GC= `grep '%GC' ./$file/fastqc_data.txt`;
$GC=(split(/\s+/,$GC))[1];
chomp $GC;

# 提取Q20，Q30
## 读入Per sequence quality scores部分的信息，保存成哈希
open FH , "<./$file/fastqc_data.txt";
while (<FH>)
    {
    next unless /#Quality/;
    while (<FH>)
        {
        @F=split;
        $hash{$F[0]}=$F[1];
        last if />>END_MODULE/;
        }
    }
## 统计Q20，Q30
$all=0;$Q20=0;$Q30=0;
$all+=$hash{$_} foreach keys %hash;
$Q20+=$hash{$_} foreach 0..20;
$Q30+=$hash{$_} foreach 0..30;
$Q20=1-$Q20/$all;
$Q30=1-$Q30/$all;
print "$file\t$total_reads\t$GC\t$Q20\t$Q30\n";
}

执行方法：

1
$ perl fastqc_stat.pl

Python

1. 格式化FASTA文件

在规范的Fasta文件中，你看到的一条序列记录包括两部分：

  • 以>起始的序列名称，占一行；
  • 由核苷酸ATCGN（核酸序列）或氨基酸字符组成的字符串，一般每60个字符一行，若一条序列很长，那么它可能会占多行；

有的时候因为一些原因（一般都是自己在上游分析时图方便生成的）你得到的fasta文件中的序列部分没有按照60个字符一行的形式进行组织，而是将整条序列放在一行里，虽然一般来说这并不会对你的分析产生太大的影响，但是进行查看的时候会有一些不方便，比如

这个时候如果想将它调整组成规范的格式，要怎么实现呢？

用BioPython将原始FASTA文件读入，然后在写出，就能得到你想要的效果

这个脚本很简单，只有四行代码

脚本名：formatFasta.py

1
2
3
4
5
6
from Bio import SeqIO
import sys

Seq = [ seq for seq in SeqIO.parse(sys.argv[1],'fasta')]

SeqIO.write(Seq,sys.argv[2],'fasta')

用法：

1
$ python formatFasta.py <in.fa> <out.fa>

R

1. 合并多个样本的定量结果成一个大矩阵 (profile matrix)

在RNAseq或者metagenome shotgun中，需要对每个样本进行逐一定量，生成每个样本各种的定量结果，然后需要合并每个样本的定量结果，形成一个包含所有样本定量的profile matrix

我们这个脚本的要实现的就是合并多个样本的定量结果成一个大矩阵 (profile matrix)

实现思路：

  • 想要保证所有样本的定量结果文件都保存在一个文件夹下，且每个文件的命名形式为：sample+suffix，即样本名+固定后缀；

  • 每个文件中只有两列信息，第一列为feature（对于RNAseq，它的feature是gene或transcript；对于metagenome shotgun，它的feature是各个物种分类层级），第二列是定量结果（quant，对于RNAseq，它的定量可以是reads count、FPKM/RPKM、TPM等；对于metagenome shotgun，它的定量一般是一个取值在[0,1]的比例值）；

  • 根据指定的文件夹（下面记作dir）和文件后缀（下面记作pattern），将文件路径为dir/*pattern的文件逐一读入，获取它们的第一列，即features的列表，从而得到total unique features的列表，同时记下样本列表（样本名为：文件名-pattern）；

  • 根据上一步得到的total unique features列表和样本列表，创建一个行数等于total unique features列表长度，列数等于样本列表的矩阵（记作$\text{MergeMatrix}_{ij}$，i表示$\text{feature}_i$，j表示$\text{sample}^j$），行名对应total unique features列表，列名对应样本列表，矩阵中每个元素的值初始化为0；

  • 再一次将文件路径为dir/*pattern的文件逐一读入，更新MergeMatrix对应列的值：若当前读入的文件名为sample_j.pattern，则当前样本为$\text{sample}^j$，则根据$\text{MergeMatrix}_{ij}=\text{sample}i^j$，对矩阵$\text{MergeMatrix}{ij}$进行更新；

脚本名：Merge2ProfileMatrix.R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
loadMatrix <- function(dir,pattern){
	files <- list.files(dir)
	clones <- vector()	# 用于保存所有样本中出现的克隆
	samples <- vector()	# 用于保存所有样本的id
	# 获取所有样本中出现的unique克隆，与所有样本的id
	for (file in files) {
		if (grepl(pattern,file)){
			data <- read.table(paste(dir,file,sep="/"),header=F,sep="\t")
			clones <- c(clones,as.character(data$V1))
			samples <- c(samples,sub(pattern,"",file))
		}
	}
	# 创建用于保存所有样本整合的表达谱，以0填充
	clones <- unique(clones)
	dataMat <- matrix(rep(0,length(samples)*length(clones)),nrow = length(clones),ncol = length(samples))
	rownames(dataMat) <- clones
	colnames(dataMat) <- samples
	# 逐一读入样本的count文件，对matrix中对应的元素进行赋值
	for (file in files) {
		if (grepl(pattern,file)){
			data <- read.table(paste(dir,file,sep="/"),header=F,sep="\t")
			sample <- as.character(sub(pattern,"",file))
			index <- match(data$V1,rownames(dataMat))
			dataMat[index,sample] <- data$V2
		}
	}
	dataMat
}

该脚本只定义了一个函数，可以在R交互模式下，载人这个脚本中的函数，然后调用它，或者直接对这个脚本进行修改，然后用Rscript直接执行修改后的脚本

2. 对profiles进行wilcox检验（并行化）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# 四个参数，按顺序分别为：
# - profiles matrix file
# - sample group file 
# - outfile name
# - threads number（默认为服务器总线程-10）

Args <- commandArgs(T)

library(parallel)

# Initiate cluster
if(is.na(Args[4])) {
	no_cores <- detectCores() - 10
} else {
	no_cores <- as.integer(Args[4])
}
cl <- makeCluster(no_cores)

data <- read.table(Args[1],head=T,row.names=1,sep='\t')
group <- read.table(Args[2],head=T,sep='\t')

# 找出profiles中每列代表的样本所属的组
col_group1 <- colnames(data) %in% sub('-','.',paste(group$Sample[group$Group==1],'.rpkm',sep=''))
col_group2 <- colnames(data) %in% sub('-','.',paste(group$Sample[group$Group==2],'.rpkm',sep=''))

# 对单个基因进行wilcox检验，并输出统计检验结果，输出格式如下：
# 	geneName,pvalue,group1Median,group2Median,direction
wilcox_fun <- function(data,col_group1,col_group2){
	g1 <- as.numeric(data[col_group1])
	g2 <- as.numeric(data[col_group2])
	# 执行wilcox检验
	stat <- wilcox.test(g1,g2,paired = FALSE, exact=NULL, correct=TRUE, alternative="two.sided")
	if(is.na(stat$p.value)){
		stat$p.value <- 1
	}
	# 对有统计学意义的基因进行判断，是上调"up"还是下调"down"或者是不变"-"（对于组2，即不吃药组）
	if(stat$p.value < 0.1  & median(g1) < median(g2)){
		G12 <- c(ifelse(is.na(stat$p.value),1,stat$p.value),median(g1),median(g2),'down')
	}
	else if(stat$p.value < 0.1  & median(g1) > median(g2)){
		G12 <- c(ifelse(is.na(stat$p.value),1,stat$p.value),median(g1),median(g2),'up')
	}
	else if(stat$p.value < 0.1  & median(g1) == median(g2)){
		G12 <- c(ifelse(is.na(stat$p.value),1,stat$p.value),median(g1),median(g2),'-')
	}
	else{
		G12 <- c()
	}
	G12
}
	
statOut <- parApply(cl,data,1,wilcox_fun,col_group1,col_group2)
stopCluster(cl)
# 删除为NULL的列表元素
for(i in names(statOut)){
	if(is.null(statOut[[i]])){
		statOut[[i]] <- NULL
	}
}

# 将结果写入文件中
if(!is.null(statOut)){
	# 转换成数据框
	statOut <- as.data.frame(t(as.matrix(as.data.frame(statOut))))
	colnames(statOut) <- c('pvalue','group1Median','group2Median','direction')
	# 写入文件
	write.table(statOut,Args[3],sep = '\t',row.names = T,col.names = T,quote = F)
}

3. 对差异显著的subprofile进行可视化——箱线图与热图

想法：

从上游的分析中得到两种样本profile的差异统计检验结果，其中有p-value和p-adjust，可以根据p-value或p-adjust设置阈值，筛选出有统计学显著性的差异observation

根据这些observation的Id从原始profile中将这些差异observation的部分提取了得到subprofile

对subprofile进行可视化——箱线图（boxplot）与热图（heatmap）

需要提供的输入：

  • 上游的分析中得到两种样本profile的差异统计检验结果

  • 原始profile

  • 样本的分组信息

脚本名：plotDiffSubprofile.R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
# 实现功能说明：导入差异分析的结果与原始profile，提取原始profile中提取
#	差异显著（通过第3个参数设置筛选的阈值）的集合绘制热图/箱线图
#
#	考虑到可能直接得到的画图输出不满足用户的需求，所以提供了一个参数来控
#	制是否输出画图用的数据，以便用户可以直接载人绘图数据，根据自己的需求
#	灵活画图
#########################################################################
# 必须参数说明：
# - output file from wilcox stat
# - profiles matrix file
# - padj threshold
# - min memmbers for subprofile
# - sample group file, 2 columns, 1st column corresponding to sample, 2nd
#	 column corresponding to group
# - logic value, False or True, whether to output prepared data frame for
#	 ggplot-boxplot and filted subprofile for heatmap, saved as raw txt and
#	 Rdata(subprofile_melt--boxplot, subprofile_matrix,anno_col--heatmap) formats
#########################################################################

library(reshape2)
library(ggplot2)
library(pheatmap)

Args <- commandArgs(T)

stat <- read.table(Args[1],header=T,sep='\t')
profile <- read.table(Args[2],head=T,sep='\t')
padj <- as.numeric(Args[3])
num_min <- as.numeric(Args[4])
group <- read.table(Args[5],head=T,sep='\t')

stat_filt <- stat[stat$pAdj<=padj,]

# 若差异的集合太小则不进行后续的绘图操作
if(dim(stat_filt)[1]>=num_min){
	# 提取出差异的subprofile
	subprofile <- profile[profile[,1] %in% rownames(stat_filt),]
	
	# 将数据框的短格式展开为长格式
	subprofile_melt <- melt(subprofile,id.vars=1,variable.name = 'sample',value.name = 'abundance')
	colnames(subprofile_melt) <- c('ID','sample','abundance')
	
	# 添加分组信息
	subprofile_melt$group <- NULL
	group_factor <- unique(group[,2])	# 获得所有组别及其表示符
	for(i in group_factor){
		subprofile_melt$group[subprofile_melt$sample %in% sub('-','.',group[group[,2]==i,1])] <- i
	}
	subprofile_melt$group <- as.factor(subprofile_melt$group) # 分组变量必须为factor
	
	# 画boxplot
	# png(paste(Args[2],".boxplotDiff.png",sep=''))
	p_box <- ggplot(subprofile_melt)
	p_box + geom_boxplot(aes(x=ID,y=abundance,fill=group))
	# dev.off()
	ggsave(paste(Args[2],".boxplotDiff.png",sep=''))
	
	# 画heatmap
	## 整理绘图数据
	subprofile_matrix <- as.matrix(subprofile[,-1])
	rownames(subprofile_matrix) <- subprofile[,1]
	## 准备heatmap的注释信息
	anno_col <- data.frame(group=as.factor(group[,2]))
	rownames(anno_col) <- sub('-','.',group[,1])
	png(paste(Args[2],".heatmapDiff.png",sep=''))
	pheatmap(subprofile_matrix,scale="row",cluster_rows=T,cluster_cols=T,annotation_col=anno_col,show_rownames = F)
	dev.off()
	
	# 若最后一个参数为True则写出数据
	if(Args[6]){
		write.table(subprofile_melt,paste(Args[2],".boxplotDiff.data",sep=''),row.names=F,col.names=T,sep='\t',quote=F) # 画箱线图的数据
		write.table(subprofile_matrix,paste(Args[2],".heatmapDiff.data",sep=''),row.names=T,col.names=T,sep='\t',quote=F) # 画热图的数据
		write.table(anno_col,paste(Args[2],".heatmapDiff.anno",sep=''),row.names=T,col.names=T,sep='\t',quote=F) # 画热图的注释数据
		save(subprofile_melt,subprofile_matrix,anno_col,file=paste(Args[2],".Rdata",sep=''))
	}
}else{
	print("Not enough memmbers of subprofile passing the padj threshold!")
}

4. 对大矩阵计算Jaccard Index（并行化+分治）

在对一个大矩阵执行相关性计算或Jaccard Index的计算时，其实执行的是矩阵任意两行（这里假设要进行分析的对象是矩阵的每个行）之间的两两的计算，若这个矩阵的规模非常庞大，有n行时，计算的时间复杂度就是$O(n^2)$，这个时候可以采用并行化策略来加速这个进程（参考上文的 2. R中的并行化方法）：

1
StatOut <- parApply(cl, data, 1, fun, data)

这样就会实现将一个 n vs. n 的问题拆分成 n 个可以并行解决的 1 vs. n 的问题，当然通过设置线程数为$m,\,(m\le n)$，使得每次只并行执行m个 1 vs. n 的问题

然后再在函数fun内部再定义一个并行化计算来进一步并行化加速上面产生的 1 vs. n 的计算：

1
2
3
4
fun <- function(vec, data){
	...
	parApply(cl, data, 1, fun2, vec)
}

在这个函数内部实现了将一个 1 vs. n 的问题拆分成 n 个可以并行解决的 1 vs. 1 的问题

这样就实现了两步并行化，这样在保证硬件条件满足的情况下，的确能显著加快分析速度

但是并行化技术会带来一个问题是，虽然时间开销减少了，但是空间开销显著增加了

比如，第一次并行化，将一个 n vs. n 的问题拆分成 $\frac{n}{m}$ 次可以并行解决的 m个 1 vs. n 的问题，则需要在每一次并行化任务中拷贝一个 1 vs. n 的计算对象，即原始有n行的矩阵被拷贝了m次，则相应的缓存空间也增加了m倍，很显然内存的占用大大增加了

空间开销显著增加带来的后果就是，很容易导致运行内存不足程序运行中断的问题，那该怎么解决这个问题呢？

可以采用分治方法（Divide-Conquer)，将原始大矩阵，按照行拆分成多个小的子块，对每个子块执行计算，从而得到每个子块的运算结果，最后再讲每个子块的结果进行合并：

脚本名：JaccardIndex.R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

###########################################################
# 	该脚本用于对Observed矩阵（行为Feature，列为Sample）
# 对feature（即行）计算Jaccard Index，可以进行并行化计算
# 	考虑到原始矩阵可能过大（主要是行太多），如果直接对
# 它执行并行化计算的话，非常容易导致内存溢出，为了解决
# 这个问题，引入了分治方法：
# 	将原始矩阵按照行等分n份（平均每份有abs(N/n)，最后一
# 份可能无法进行等分，最后一份为N%n），然后对这n份的每两
# 份之间分别执行并行化计算，则总共有n^2次，得到n^2个子结
# 果进行合并
###########################################################

# 参数说明：
# （1）Observed Matrix，每行表示一个feature，每列表示一个样本
# （2）并行化使用的线程数
# （3）切分的分数（默认不切块）


library(parallel)


Args <- commandArgs(TRUE)

inFile <- Args[1]

# 开启并行化
if(!is.na(Args[2])){
	no_cores <- as.integer(Args[2])
	cl <- makeCluster(no_cores)
}

data <- read.table(inFile,header=T,row.names=1)

print("Load Observed Matrix Successfully")

# 计算Jaccard Index
JaccardIndexSer <- function(SourceVec, TargetMatrix){
	JaccardIndexOne <- function(SourceVec, TargetVec){
		sum(SourceVec == TargetVec)/length(SourceVec)
	}
	apply(TargetMatrix, 1, JaccardIndexOne, SourceVec)
}

# 执行分治方法
if(!is.na(Args[3])){
	n_row <- nrow(data)
	nblock <- Args[3]
	nrow_block <- ceiling(n_row/nblock)
	StatOutList <- vector("list", nblock)
	
	# 1. 开始进行分块计算
	print("[Divide-Conquer]Start carry out Divide-Conquer for Large Observed Matrix")
	print("##################################################")
	for(i in 1:nblock){
		for(j in 1:nblock){
			nrow_start <- (i-1)*nrow_block+1
			nrow_end <- i*nrow_block
			# 并行化计算
			if(!is.na(Args[2])){
				print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index parallel for block: ",i,"-",j,sep=''))
				StatOutList[[i]] <- append(StatOutList[[i]], parApply(cl, data[nrow_start:nrow_end,], 1 , JaccardIndexSer, data))
				print(paste("[Divide-Conquer]Finish run parallel for block: ",i,"-",j,sep=''))
			# 串行计算
			}else{
				print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index serially for block: ",i,"-",j,sep=''))
				StatOutList[[i]] <- append(StatOutList[[i]], apply(data, 1 , JaccardIndexSer, data))
				print(paste("[Divide-Conquer]Finish run serially for block: ",i,"-",j,sep=''))
			}
		}
	}
	# 2. 结束分治方法的分块计算
	if(!is.na(Args[2])){
		print("##################################################")
		print("[Divide-Conquer]Finish parallel running for statistic Jaccard Index!")
		stopCluster(cl)
	}else{
		print("##################################################")
		print("[Divide-Conquer]Finish serial running for statistic Jaccard Index!")
	}
	# 3. 开始进行子块结果的合并
	print("[Divide-Conquer]Start bind sub-block statout")
	StatOut <- vector("list", nblock)
	# 先对列进行合并
	for(i in 1:nblock){
		for(block in StatOutList[[i]]){
			StatOut[[i]] <- cbind(StatOut[[i]], block)
		}
	}
	# 再对行进行合并
	StatOutMerge <- data.frame()
	for(block in StatOut){
		StatOutMerge <- rbind(StatOutMerge, block)
	}
	StatOut <- StatOutMerge
# 不执行分治方法
}else{	
# 并行化计算
	if(!is.na(Args[2])){
		print("[Common]Start carry out statistic Jaccard Index parallel")
		StatOut <- parApply(cl, data, 1 , JaccardIndexSer, data)
		print("[Common]Finish run parallel")
		stopCluster(cl)
	# 串行计算
	}else{
		print("[Common]Start carry out statistic Jaccard Index serially")
		StatOut <- apply(data, 1 , JaccardIndexSer, data)
		print("[Common]Finish run serially")
	}
}
save(StatOut,'JaccardIndex.Rdata')

5. 合并两个矩阵

在什么情节下，我们需要合并两个矩阵？

对于两个批次的数据，我们分别用上面 1. 合并多个样本的定量结果成一个大矩阵 (profile matrix)得到了两个profile matrix，这两个matrix的行表示feature，列表示sample，样本一般不会有重复，而features大部分是重叠的，也有少部分是不重叠的，那么将这两个矩阵合并之后，这个新矩阵的行数为两个矩阵features的并集的大小，列数为两个矩阵列数的和

脚本名：MergeTwoMat.R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
Args <- commandArgs(T)

Matrix1 <- read.table(Args[1], header=T, row.names=1)
Matrix2 <- read.table(Args[2], header=T, row.names=1)

# 获取两个矩阵feature的并集
MergeFeatureList <- as.character(union(rownames(Matrix1),rownames(Matrix2)))
# 获取两个矩阵的sample的并集（默认没有重叠）
ColList <- as.character(c(colnames(Matrix1), colnames(Matrix2)))

# 初始化合并后的新矩阵，行数为feature并集大小，列数为sample并集大小
MergeMatrix <- matrix(rep(0,length(MergeFeatureList)*length(ColList)), nrow = length(MergeFeatureList), ncol = length(ColList))
rownames(MergeMatrix) <- MergeCloneList
colnames(MergeMatrix) <- ColList

# 对第一个矩阵中对应的列（即样本）进行更新
index <- match(rownames(Matrix1), MergeCloneList)
for(col in colnames(Matrix1)){
	MergeMatrix[index, col] <- Matrix1[,col]
}

# 对第二个矩阵中对应的列（即样本）进行更新
index <- match(rownames(Matrix2), MergeCloneList)
for(col in colnames(Matrix2)){
	MergeMatrix[index, col] <- Matrix2[,col]
}

用法：

1
$ Rscript MergeTwoMat.R <matrix1> <matrix2>