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

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

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

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

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

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

## ****************************************************************
## ** 该脚本用于将输入的表格或矩阵文件，按照列平均拆分为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];
}