Perl
Python
- 1. 格式化FASTA文件
R

实用小脚本

Perl ^目录

1. 拆分FASTA文件 ^目录

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

脚本名：splitFasta.pl

#!/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;

执行方法：

$ 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

#!/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;

执行方法：

$ 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

#!/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;

执行方法：

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

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

4. 根据序列Id提取FASTA序列 ^目录

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

脚本名：extractSeqFromFasta.pl

#!/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(/^>(.*?)\s/){
		$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 @S = split /\t/;
	my $gene = $S[0];
	if($Hash_fasta{$gene}){
		print OUT ">$gene\n$Hash_fasta{$gene}\n";
	}
}
close(LIST);
close(OUT);

使用方法：

$ 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

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";
}

执行方法：

$ perl fastqc_stat.pl

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

该功能的另一种实现方式见R语言部分的 1. 合并多个样本的定量结果成一个大矩阵 (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的文件逐一读入，然后保存为双重哈希形式，即

$$%H={feature\to{sample\to count}}$$

同时记录下样本名的列表$@S$

然后，根据上面构造出来的双重哈希$%H$和样本名列表$@S$，对每个feature逐行写出到输出文件中，若当前feature为$i$，遍历样本名列表，若当前样本名为$j$，则feature $i$ 在样本$i$的取值记为$n_{ij}$

$$n_{ij}=\quad \text{defined}(H{j})\quad ?\quad H{j}:0$$

脚本名：matrix_maker.pl

#!/usr/bin/perl -w
use strict;
use POSIX;
use Data::Dumper;
use String::Util qw(trim);
use Term::ANSIColor;

## *********************************************
## ** version history **
## *********************************************
my $ver = '1.0';
my $last_modified = 'Sep 26, 2018';

## *********************************************
## ** GET opts **
## *********************************************
use Getopt::Long;
my %opts=();
GetOptions(\%opts,"i:s","o:s", "samples:n", "reads:n", "debug");

if (!$opts{i} or !$opts{o} ){
	print "--------------------------------------------------------------------------------------------------
	\t\tversion : $ver ; last modified : $last_modified
--------------------------------------------------------------------------------------------------
	USAGE: perl $0
		-i input dir
		-o output file
	[optional]
		-samples number of samples in which a CDR3 gene should be found,
			otherwise this gene will be removed
			default is 0: does not do this removing
			logic: >= 
		-reads number of reads mapped to a CDR3 gene
			otherwise this gene will be removed (set to 0)
			default 0: does not do filtering
			logic: >= 
		-debug debug mode, if present, will only process the first 10 files ...
--------------------------------------------------------------------------------------------------\n";
	exit;
}

## --- user input parameters  ---
my $min         = defined $opts{samples} ? $opts{samples} : 0;
my $reads       = defined $opts{reads} ? $opts{reads} : 0;
my $debugmode   = defined $opts{debug};

if( $debugmode ){
	print STDERR "==========================================================\n";
	print STDERR "\tDEBUG MODE: only the first 10 files will be processed!!!\n";
	print STDERR "==========================================================\n";
}

if( $min or $reads ){
	print STDERR "====================CUTOFFs SET ==========================\n";
	print STDERR "\tmin   reads: $reads \n" if( $reads );
	print STDERR "\tmin samples: $min \n" if( $min );
	print STDERR "==========================================================\n";
}

## *********************************************
## ** Load files **
## *********************************************
my %hAA2File2Count = (); ### $hash{ $AA }{ $File } = $count;


### read dir and process --
my $file_count = 0;

print STDERR "\tread dir and load frequencies ... \n";
my @aFiles = (); 
my $indir = $opts{i}; chop $indir if ($indir =~ /\/$/);
opendir(DIR, $indir) || die "can't opendir $indir: $!";
# 只对指点输入文件夹的非隐藏文件进行处理
foreach my $infile (grep { !/^\./ && -f "$indir/$_" } readdir(DIR)){
	## keep tracking infile --
	push @aFiles, $infile;
	
	### ---- open INFILE ----
	my $infile_fullpath = $indir . '/' . $infile;
	open INFILE, $infile_fullpath or die;
	# 用双重哈希来保存读入的数据：{aa -> {infile/sample -> count}}
	while(<INFILE>){
		chomp;
		my ( $aa, $count ) = split(/\t/, $_);
		$hAA2File2Count{ $aa }{ $infile } = $count if( $count >= $reads ); # 卡最低reads数的阈值
	}
	close INFILE;
	
	$file_count ++;
	print STDERR "\t\t", $file_count, ": ", $infile , " loaded ... \n";
	
	## -- if debug mode --
	if( $debugmode and $file_count >= 10 ){
		print STDERR "\tbreak from loop due to DEBUG MODE, the results MIGHT BE IMCOMPLETE!!!\n";
		last;
	}
}
closedir DIR;
print STDERR "\tin total $file_count files loaded; \n\n";


### =============== now prepare output ======================
print STDERR "\tnow ready to generate output ... \n";
if( $min or $reads ){
	print STDERR "\tNOTE: only genes that are found in >= $min samples and supported by >= $reads reads will be included in the output file ... \n";    
}

## --- open OUT file for output --
open OUT, ">$opts{o}" or die;
print OUT join("\t", "GENE", @aFiles), "\n"; ## print the title line

## --- keep tracking genes processed ... ---
my $genes_count = 0;
my $genes_total = 0;
while( my ( $aa, $hashref ) = each %hAA2File2Count ){
	my $files_count = scalar keys %{ $hashref };
	$genes_total ++;
	
	## --- 卡最低表达样本数的阈值，来决定是否保留该基因/克隆 --- 
	if( $files_count >= $min ){
		$genes_count ++;
		my @aRow = ();
		foreach my $infile( @aFiles ){
			push @aRow, exists $$hashref{ $infile } ? $$hashref{ $infile } : 0;
		}
		
		print OUT join("\t", $aa, @aRow ), "\n";
		
		## --- print some debug information ... 
		if( $genes_count % 1000000 == 0 ){
			print STDERR "\t\t$genes_count genes written to out file, out of in total $genes_total genes ... \n";
		}
	}
}
close OUT;
print STDERR "\tin total ", &commify($genes_count), " written to output ... out of in total ", &commify($genes_total), " genes processed.\n";
print STDERR "\tin total ", sprintf("%.3f", $genes_count / $genes_total * 100), "% of genes retained!!\n";
print STDERR "\tall jobs done; \n\n";

sub commify {
	my $text = reverse $_[0];
	$text =~ s/(\d\d\d)(?=\d)(?!\d*\.)/$1,/g;
	return scalar reverse $text;
}

7. 将大矩阵文件按照固定列数进行分割 ^目录

该脚本的难度在于，通过指定固定的子矩阵文件的列数来进行原始矩阵切割，最终切割得到的子文件数是不确定的，需要根据给定的输入文件的总列数来推算出来

有两种解决方案：

多次扫描输入文件，一次产生一个子文件

扫描一次输入文件，得到一个子矩阵文件，再重新从头到尾扫描一次，得到另一个子矩阵文件……若最终得到的子文件的数量为n，则需要进行n次从头到尾的扫描——当输入文件很大，尤其是行数很多时，效率极低；
只扫描输入文件一次，同时产生多个子文件
只创建一个文件句柄，让这个句柄不断地开启和关闭，一次开启和关闭算一轮，每从原文件中读入一行，然后分列写入不同子文件时，需要每一轮让其指向不同的文件，核心代码如下：
 while(<I>){
 	chomp;
 	@row = split /\t/;
 	$feature_name = shift @row;
 	for($i=0;$i<$n_subfiles;$i++){
 		open O,">>${outpreffix}_$i" or die "Can not create file: ${outpreffix}_$i!\n";
 		print O "$feature_name\t";
 		print O join("\t",@row[($i*$n_col)..(($i + 1)*$n_col -1)])."\n";
 		close O;
 	}
 }
该方案的问题在于，需要不断地开关文件句柄，而文件句柄的开关操作是有比较高的计算复杂度的——子文件越多，则文件句柄的开关轮次越多，则效率越低；
创建多个文件句柄（句柄数量和目标子文件数量相同），一个文件句柄对应一个目标子文件，则每从原文件中读入一行，然后分列写入不同子文件时，只需要不断切换其对应的文件句柄即可，不需要像上面的方法一样，进行文件句柄的开关操作

这种方法的问题在于：传统的文件句柄的创建方式为 open HANDLE,">$file" 这样的语法结构，其中的HANDLE就是文件句柄，而且是裸字文件句柄，在这种用法下，句柄名是要以确定的字符串的形式指定的，而我们要创建的子文件的数量不定，不可能进行确定的裸字文件句柄的指定

这个时候，就可以考虑使用文件句柄的另一种创建方式——变量文件句柄，即 open $fh,">$file"，然后将每一次创建得到的文件句柄追加保存到一个数组中，即 push(@FH,$fh)

不过需要注意，若以这样的方式循环创建多个文件句柄，而每次循环改变的只是最后的$file变量对应的文件路径，那么每一次得到的文件句柄会覆盖上一个文件句柄，因此虽然最后得到了保存在数组中的多个文件句柄，但是它们都指向最后的那个文件——将前面的命令改成open my $fh,">$file"就能完美地解决上面的问题，即在变量文件句柄$fh之前加上my，申明它是一个局部变量，则新创建的文件句柄虽然于上一个文件句柄拥有相同的变量名，但是它们在保存在不同的内存空间里，所以本质上不是一个变量，也就不会产生前面的新变量覆盖就变量的问题了

脚本名：SplitByCol.py

## ****************************************************************
## ** 该脚本用于将输入的表格或矩阵文件，按照列平均拆分为n个子文件
## ****************************************************************


## （1）原始的profile matrix文件
## （2）指定生成的子文件的数量
## （3）输出文件的前缀，输出文件的文件名形式为：*_1, *_2, ...
## （4）每处理n行输出一次进度报告

use POSIX;
use Getopt::Long;

# 脚本帮助文档
=head1 Description

Thise script is used to split large matrix file into serveral sub-matrix files by column

=head1 Usage

$0 -i <input> [-N <subfiles number> | -n <col number per file>] -v <row number>

=head1 Parameters

-i	[str]	Input raw matrix file
-N	[int]	Specify the number of output sub-matrix files(default: 10)
-n	[int]	Specify the column number per sub-matrix file, attention: not set "-N" and "-n" at the same time, otherwise, it will stop runing and throw out an error
-o	[str]	Output file preffix
-v	[int]	The number of processed rows per report

=cut

my ($infile, $n_subfiles, $ncol_perfile, $outpreffix, $n_report);

GetOptions(
"i:s"=>\$infile,
"N:s"=>\$n_subfiles,
"n:i"=>\$ncol_perfile,
"o:s"=>\$outpreffix,
"v:i"=>\$n_report
);

die `pod2text $0` if (!defined($infile)) or (!defined($outpreffix));

# 若同时设置"-N"和"-n"参数，则报错并停止运行
if(defined($n_subfiles) && defined($ncol_perfile)){
	print STDERR "Not set \"-N\" and \"-n\" at the same time\n\n";
	die `pod2text $0`;
}

# 若"-N"和"-n"参数都未设置，则设置默认的"-N"参数
unless(defined($n_subfiles) && defined($ncol_perfile)){
	$n_subfiles = 10;
}

# 计算对应的子文件数或子文件的列数
open I,"<$infile" or die "No such file: $infile!\n";
$header = <I>;
chomp $header;
@Header = split /\t/,$header;
$rowname = shift @Header;
$total_col = $#Header - 1; # 计算原始文件除第一列之外的总列数
# 若设置了"-N"参数，则计算每个子文件的列数
if(defined $n_subfiles){
	$ncol_perfile = ceil($total_col / $n_subfiles); # 计算子文件的平均列数

# 若设置了"-n"参数，则计算子文件的数量
}else{
	$n_subfiles = ceil($total_col / $ncol_perfile); # 计算子文件的数量
}

# 新建若干个输出文件的文件句柄，并追加保存到一个数组变量中
for($i=0;$i<$n_subfiles;$i++){
	open my $fh,">${outpreffix}_$i" or die "Can not create file: ${outpreffix}_$i!\n";
	# 写入文件的表头
	print {$fh} "$rowname\t";
	print {$fh} join("\t",@Header[($i*$ncol_perfile)..(($i + 1)*$ncol_perfile -1)])."\n";
	push(@FH,$fh);
}

# 写入除表头外的具体内容
$j = 0;
while(<I>){
	chomp;
	# 每处理1万行，输出一次进度报告
	$j++;
	if($j%$n_report==0){
		print STDERR "$j line processed\n";
	}
	@row = split /\t/;
	$feature_name = shift @row;
	for($i=0;$i<$n_subfiles;$i++){
		print {$FH[$i]} "$feature_name\t";
		print {$FH[$i]} join("\t",@row[($i*$ncol_perfile)..(($i + 1)*$ncol_perfile -1)])."\n";
	}
}

print STDERR "Finish spliting raw matrix file\n";
for($i=0;$i<$n_subfiles;$i++){
	close $FH[$i];
}

用法：

perl SplitByCol.py -i <input> [-N <subfiles number> | -n <col number per file>] -v <row number>

Python ^目录

1. 格式化FASTA文件 ^目录

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

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

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

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

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

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

脚本名：formatFasta.py

from Bio import SeqIO
import sys

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

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

用法：

$ 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

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检验 ^目录

# 四个参数，按顺序分别为：
# - 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

########################################################################
# 实现功能说明：导入差异分析的结果与原始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中的并行化方法）：

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 的计算：

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


###########################################################
# 	该脚本用于对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

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]
}

用法：

$ Rscript MergeTwoMat.R <matrix1> <matrix2>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

实用小脚本.md

实用小脚本.md

实用小脚本

Perl ^目录

1. 拆分FASTA文件 ^目录

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

3. 双端FASTQ文件进行双端配对 ^目录

4. 根据序列Id提取FASTA序列 ^目录

5. 统计FastQC输出 ^目录

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

7. 将大矩阵文件按照固定列数进行分割 ^目录

Python ^目录

1. 格式化FASTA文件 ^目录

R ^目录

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

2. 对profiles进行wilcox检验 ^目录

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

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

5. 合并两个矩阵 ^目录

Files

实用小脚本.md

Latest commit

History

实用小脚本.md

File metadata and controls

实用小脚本

Perl 目录

1. 拆分FASTA文件 目录

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

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

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

5. 统计FastQC输出 目录

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

7. 将大矩阵文件按照固定列数进行分割 目录

Python 目录

1. 格式化FASTA文件 目录

R 目录

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

2. 对profiles进行wilcox检验 目录

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

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

5. 合并两个矩阵 目录

Perl ^目录

1. 拆分FASTA文件 ^目录

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

3. 双端FASTQ文件进行双端配对 ^目录

4. 根据序列Id提取FASTA序列 ^目录

5. 统计FastQC输出 ^目录

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

7. 将大矩阵文件按照固定列数进行分割 ^目录

Python ^目录

1. 格式化FASTA文件 ^目录

R ^目录

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

2. 对profiles进行wilcox检验 ^目录

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

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

5. 合并两个矩阵 ^目录