Arichard:Computational/fastq2CNV.pl
- !/usr/bin/perl -w
use strict;
my $system; my $project_info_file = $ARGV[0]; my $aln_core = 8; open(INFILE, "$project_info_file")||die("Error in opening file $project_info_file\n");
- Read in system
- my $line = <INFILE>;
- chop($line);
- $system = $line;
$system = "genome-miner";
- Read Header Line
my $line = <INFILE>;
- System-specific parameters
my ($mapper, $samtool, $picardDir, $SNS_varbin_counter, $bowtie_template, $samtool_template); if ($system eq "genome-miner") { $aln_core = 4; $mapper = "/home/kunzhang/softwares/bowtie-latest/bowtie"; $samtool = "/home/kunzhang/softwares/samtools-latest/samtools"; $picardDir="java -Xmx4g -jar /home/kunzhang/softwares/picard-tools-latest/"; $bowtie_template = "/GenomeDB/UCSC/Homo_sapiens/hg19/Sequence/BowtieIndex/genome"; $samtool_template = "/GenomeDB/UCSC/Homo_sapiens/hg19/Sequence/BowtieIndex/genome.fa"; $SNS_varbin_counter = "/home/kunzhang/softwares/CSHL_SNS/varbin.50k.sam.triton.py"; } elsif ($system eq "triton") { $mapper = "/home/k4zhang/softwares/bowtie-latest/bowtie"; $samtool = "/projects/zhang-lab/softwares/samtools-0.1.18/samtools"; $picardDir="java -Xmx4g -jar /projects/zhang-lab/softwares/picard-tools-latest/"; $bowtie_template = "/home/k4zhang/LabDir/RNAseq/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex/genome"; $samtool_template = "/home/k4zhang/LabDir/RNAseq/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa"; $SNS_varbin_counter = "/home/k4zhang/softwares/CSHL_SNS/varbin.50k.sam.triton.py"; } else { die("Please choose a system (genome-miner or triton) on line 1 of the input file.\n"); }
my @words = split(/\./, $project_info_file); my $project_handle=$words[0]; my %sampleTable; while($line = <INFILE>){ chomp($line); my @fields = split(/[\t ]/, $line); next if(scalar(@fields)<5); my $key = $fields[2] . "_" . $fields[3] . "_" . $fields[4]; $sampleTable{$key}->{"SM"}=$fields[0]; $sampleTable{$key}->{"LIB"}=$fields[1];
if (!(-e $fields[5])) { die("Could not find file $fields[5]\n"); } $sampleTable{$key}->{"Fq1"}=$fields[5];
} close(INFILE);
my $cmd; my @args;
my $fastq_se = $project_handle . ".1.fq";
my @merge_bam_files; foreach my $key (sort keys(%sampleTable)){ $cmd=sprintf("less %s > %s", $sampleTable{$key}->{"Fq1"}, $fastq_se); print "$cmd\n"; system($cmd); my $sam_file = $key . ".sam"; my $bam_file = $key . ".bam";
$cmd = "$mapper -p $aln_core --sam-RG \"ID:$key\" --sam-RG \"PL:ILLUMINA\" --sam-RG \"LB:".$sampleTable{$key}->{"LIB"}."\" --sam-RG \"SM:".$sampleTable{$key}->{"SM"}."\" -S -t -n2 -e 70 -m 1 --phred64-quals --best --strata $bowtie_template $fastq_se $sam_file"; print "$cmd\n"; system($cmd); $cmd ="$samtool view -b -t $samtool_template $sam_file > $bam_file"; print "$cmd\n"; system($cmd); push(@merge_bam_files, $bam_file); #unlink($sam_file); unlink($fastq_se); }
my $project_bam_file = $project_handle . ".bam";
my $joinString = " "; for(my $counter = 0;$counter<@merge_bam_files;$counter++) { $joinString = $joinString." INPUT=".$merge_bam_files[$counter]; } $cmd = $picardDir . "MergeSamFiles.jar TMP_DIR=./tmp/ $joinString OUTPUT=$project_bam_file VALIDATION_STRINGENCY=SILENT"; #CHANGED print "$cmd\n"; system($cmd); #unlink(@merge_bam_files);
system("mkdir ./tmp");
my $sorted_bam_file = $project_handle . ".sorted.bam"; $cmd = $picardDir . "SortSam.jar TMP_DIR=./tmp/ INPUT=$project_bam_file OUTPUT=$sorted_bam_file QUIET=True SORT_ORDER=coordinate VALIDATION_STRINGENCY=SILENT"; print "$cmd\n"; system($cmd);
my $unique_bam_file = $project_handle . ".unique.bam"; $cmd = $samtool . " rmdup -S $sorted_bam_file $unique_bam_file"; print "$cmd\n"; system($cmd); unlink($sorted_bam_file);
unlink($project_bam_file);
#Creating bam index
$cmd = "$samtool index $unique_bam_file";
print "$cmd\n"; system($cmd);
#Generate a sam file
- my $unique_sam_file = $project_handle . ".unique.sam";
- $cmd = "$samtool view -o $unique_sam_file $unique_bam_file";
- print "$cmd\n";
- system($cmd);
#read counting in variable bins
- my $var_bin_counting_out=$project_handle.".varbin50k.out";
- my $var_bin_counting_stat=$project_handle.".varbin50k.stat.txt";
- $cmd = "$SNS_varbin_counter $unique_sam_file $var_bin_counting_out $var_bin_counting_stat";
- print "$cmd\n";
- system($cmd);
- unlink($unique_sam_file);