Arichard:Computational/fastq2CNV.pl

From ZhangLabWiki
Jump to navigation Jump to search
  1. !/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");

  1. Read in system
  2. my $line = <INFILE>;
  3. chop($line);
  4. $system = $line;

$system = "genome-miner";

  1. Read Header Line

my $line = <INFILE>;

  1. 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

  1. my $unique_sam_file = $project_handle . ".unique.sam";
  2. $cmd = "$samtool view -o $unique_sam_file $unique_bam_file";
  3. print "$cmd\n";
  4. system($cmd);

#read counting in variable bins

  1. my $var_bin_counting_out=$project_handle.".varbin50k.out";
  2. my $var_bin_counting_stat=$project_handle.".varbin50k.stat.txt";
  3. $cmd = "$SNS_varbin_counter $unique_sam_file $var_bin_counting_out $var_bin_counting_stat";
  4. print "$cmd\n";
  5. system($cmd);
  6. unlink($unique_sam_file);