classifiers

http://www.cbs.dtu.dk/~hanni/ALL/Supplement(R%20packages%20and%20functions).htm

How to find the best K (# of clusters) in consensus clustering

It is the start point of the line, which drops the most sharply (i.e. max negative slope).
e.g. K=3 (in green circle)

Endnote Journal Style

http://endnote.com/downloads/styles

search by journal name

copy to C:\Program Files (x86)\EndNote X3\Styles

NSG Illumina pipeline screen output

Illumina 1.3+ fastq format: ASCII(min, max) = (66, 102)
2012/08/25 11:41:15 START maq ill2sanger Run1_testicular-28T_lane2_read1_sequence.txt Testis_T28_read1_sanger.fq
2012/08/25 11:42:57 SUCCESS after running 0 hours 1 minutes 42 seconds

2012/08/25 11:42:57 START maq ill2sanger Run1_testicular-28T_lane2_read2_sequence.txt Testis_T28_read2_sanger.fq
2012/08/25 11:44:37 SUCCESS after running 0 hours 1 minutes 40 seconds

2012/08/25 11:44:37 START fastx_quality_stats -Q 33 -i Testis_T28_read1_sanger.fq -o ./qcstats/Testis_T28_read1_sanger.fq_stats.txt
2012/08/25 11:49:13 SUCCESS after running 0 hours 4 minutes 36 seconds

2012/08/25 11:49:13 START fastq_quality_boxplot_graph.sh -i ./qcstats/Testis_T28_read1_sanger.fq_stats.txt -o ./qcstats/Testis_T28_read1_sanger.fq_stats.png -t "Testis_T28_read1_sanger.fq"
2012/08/25 11:49:13 SUCCESS after running 0 hours 0 minutes 0 seconds

2012/08/25 11:49:13 START fastx_quality_stats -Q 33 -i Testis_T28_read2_sanger.fq -o ./qcstats/Testis_T28_read2_sanger.fq_stats.txt
2012/08/25 11:53:51 SUCCESS after running 0 hours 4 minutes 38 seconds

2012/08/25 11:53:51 START fastq_quality_boxplot_graph.sh -i ./qcstats/Testis_T28_read2_sanger.fq_stats.txt -o ./qcstats/Testis_T28_read2_sanger.fq_stats.png -t "Testis_T28_read2_sanger.fq"
2012/08/25 11:53:51 SUCCESS after running 0 hours 0 minutes 0 seconds

2012/08/25 11:53:51 START bwa aln -t 8 /data/public/reference_genomes/genome_hg18.fa Testis_T28_read1_sanger.fq > ./bam/Testis_T28.read1.sai
2012/08/25 12:42:44 SUCCESS after running 0 hours 48 minutes 53 seconds

2012/08/25 12:42:44 START bwa aln -t 8 /data/public/reference_genomes/genome_hg18.fa Testis_T28_read2_sanger.fq > ./bam/Testis_T28.read2.sai
2012/08/25 13:33:10 SUCCESS after running 0 hours 50 minutes 26 seconds

2012/08/25 13:33:10 START bwa sampe /data/public/reference_genomes/genome_hg18.fa ./bam/Testis_T28.read1.sai ./bam/Testis_T28.read2.sai Testis_T28_read1_sanger.fq Testis_T28_read2_sanger.fq > ./bam/Testis_T28.sam
2012/08/25 13:50:34 SUCCESS after running 0 hours 17 minutes 24 seconds

2012/08/25 13:50:34 START samtools import /data/public/reference_genomes/genome_hg18.fa.fai ./bam/Testis_T28.sam ./bam/Testis_T28.bam
2012/08/25 14:03:10 SUCCESS after running 0 hours 12 minutes 36 seconds

2012/08/25 14:03:10 START samtools sort ./bam/Testis_T28.bam ./bam/Testis_T28.sorted
2012/08/25 14:25:58 SUCCESS after running 0 hours 22 minutes 48 seconds

2012/08/25 14:25:58 START samtools index ./bam/Testis_T28.sorted.bam
2012/08/25 14:26:42 SUCCESS after running 0 hours 0 minutes 44 seconds

2012/08/25 14:26:42 START java -Xmx4g -jar /home/gmslz/tools/picard-tools-1.44/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT M=PCR_duplicates.M.txt REMOVE_DUPLICATES=true AS=true I=./bam/Testis_T28.sorted.bam O=./bam/Testis_T28.sorted.rmdup.bam
2012/08/25 14:40:45 SUCCESS after running 0 hours 14 minutes 3 seconds

2012/08/25 14:40:45 START java -Xmx4g -jar /home/gmslz/tools/picard-tools-1.44/AddOrReplaceReadGroups.jar RGID=lane1 RGLB=Testis_T28.fastq RGPL=Illumina RGPU=run RGSM=Testis_T28 CREATE_INDEX=TRUE VALIDATION_STRINGENCY=SILENT I=./bam/Testis_T28.sorted.rmdup.bam O=./bam/Testis_T28.sorted.rmdup.rg.bam
2012/08/25 14:52:13 SUCCESS after running 0 hours 11 minutes 28 seconds

2012/08/25 14:52:13 START samtools index ./bam/Testis_T28.sorted.rmdup.rg.bam
2012/08/25 14:52:59 SUCCESS after running 0 hours 0 minutes 46 seconds

2012/08/25 14:53:00 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -l INFO -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.sorted.rmdup.rg.bam -o ./snp_analysis_gatk/Testis_T28.intervals > ./log/Testis_T28.intervals.log
2012/08/25 16:29:57 SUCCESS after running 1 hours 36 minutes 57 seconds

2012/08/25 16:29:57 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T IndelRealigner -l INFO -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.sorted.rmdup.rg.bam -targetIntervals ./snp_analysis_gatk/Testis_T28.intervals -o ./bam/Testis_T28.realigned.bam -stats ./bam/Testis_T28.realign.stats.txt > ./log/Testis_T28.realign.log
2012/08/26 00:40:15 SUCCESS after running 8 hours 10 minutes 18 seconds

2012/08/26 00:40:15 START samtools index ./bam/Testis_T28.realigned.bam
2012/08/26 00:41:08 SUCCESS after running 0 hours 0 minutes 53 seconds

2012/08/26 00:41:08 START java -Xmx4g -jar /home/gmslz/tools/picard-tools-1.44/SortSam.jar SO=coordinate VALIDATION_STRINGENCY=SILENT I=./bam/Testis_T28.realigned.bam O=./bam/Testis_T28.realigned.sorted.bam
2012/08/26 00:52:36 SUCCESS after running 0 hours 11 minutes 28 seconds

2012/08/26 00:52:36 START samtools index ./bam/Testis_T28.realigned.sorted.bam
2012/08/26 00:53:21 SUCCESS after running 0 hours 0 minutes 45 seconds

2012/08/26 00:53:21 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /home/gmslz/tools/picard-tools-1.44/FixMateInformation.jar SO=coordinate VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=2000000 I=./bam/Testis_T28.realigned.sorted.bam O=./bam/Testis_T28.realigned.sorted.fixed.bam
2012/08/26 01:14:28 SUCCESS after running 0 hours 21 minutes 7 seconds

2012/08/26 01:14:28 START samtools index ./bam/Testis_T28.realigned.sorted.fixed.bam
2012/08/26 01:15:13 SUCCESS after running 0 hours 0 minutes 45 seconds

2012/08/26 01:15:13 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T CountCovariates -l INFO -nt 8 --downsample_to_coverage 10000 --default_platform Illumina -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.realigned.sorted.fixed.bam -recalFile ./snp_analysis_gatk/Testis_T28.recal.csv > ./log/Testis_T28.count.covariates.log
2012/08/26 01:25:00 SUCCESS after running 0 hours 9 minutes 47 seconds

2012/08/26 01:25:00 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T TableRecalibration -l INFO --default_platform Illumina -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.realigned.sorted.fixed.bam -recalFile ./snp_analysis_gatk/Testis_T28.recal.csv --out ./bam/Testis_T28.final.recalibrated.bam > ./log/Testis_T28.recalibration.log
2012/08/26 02:00:22 SUCCESS after running 0 hours 35 minutes 22 seconds

2012/08/26 02:00:22 START samtools index ./bam/Testis_T28.final.recalibrated.bam
2012/08/26 02:01:27 SUCCESS after running 0 hours 1 minutes 5 seconds

2012/08/26 02:01:27 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T IndelGenotyperV2 -l INFO --window_size 450 -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.final.recalibrated.bam -o ./snp_analysis_gatk/Testis_T28.indels.vcf -bed ./snp_analysis_gatk/Testis_T28.indels.bed > ./log/Testis_T28.indels.log
2012/08/26 02:40:46 SUCCESS after running 0 hours 39 minutes 19 seconds

2012/08/26 02:40:46 START python /data/public/tools/gatk-v4333/Sting/python/makeIndelMask.py ./snp_analysis_gatk/Testis_T28.indels.bed 10 ./snp_analysis_gatk/Testis_T28.indels.mask.bed
2012/08/26 02:40:46 SUCCESS after running 0 hours 0 minutes 0 seconds

2012/08/26 02:40:46 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T UnifiedGenotyper -l INFO -stand_call_conf 30 -stand_emit_conf 10 -pl SOLEXA -mmq 30 -mm40 3 -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -I ./bam/Testis_T28.final.recalibrated.bam -o ./snp_analysis_gatk/Testis_T28.raw.snps.vcf > ./log/Testis_T28.snp.log
2012/08/26 04:48:58 SUCCESS after running 2 hours 8 minutes 12 seconds

2012/08/26 04:48:58 START java -Xmx4g -Djava.io.tmpdir=/home/tmp/ -jar /data/public/tools/gatk-v4333/Sting/dist/GenomeAnalysisTK.jar -T VariantFiltration -l INFO --clusterSize 3 --clusterWindowSize 10 --filterExpression "DP <= 8" --filterName "DP8" --filterExpression "SB > -0.10" --filterName "StrandBias" --filterExpression "HRun > 8" --filterName "HRun8" --filterExpression "QD < 5.0" --filterName "QD5" --filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filterName "hard2validate" -D /data/nextgen1/pipeline/dbsnp_130_hg18.rod -R /data/public/reference_genomes/genome_hg18.fa -B:variant,VCF ./snp_analysis_gatk/Testis_T28.raw.snps.vcf -B:mask,Bed ./snp_analysis_gatk/Testis_T28.indels.mask.bed --maskName InDel -o ./snp_analysis_gatk/Testis_T28.snp.filtered.vcf > ./log/Testis_T28.snp.filter.log
2012/08/26 04:53:29 SUCCESS after running 0 hours 4 minutes 31 seconds

NGS Illumina pipeline main caller

#!/usr/bin/perl
use lib './';
use NGS_Illumina;


# Assume the paths for fastx_tk, maq, bwa, samtools are in ENV{PATH}; check $PATH setting $HOME/.bash_profile
# Assume perl and python are installed
my $NGS_obj = new NGS_Illumina ('sample_name' => 'Testis_T28',
'input_fq1' => 'Run1_testicular-28T_lane2_read1_sequence.txt',
'input_fq2' => 'Run1_testicular-28T_lane2_read2_sequence.txt',
'pe' => 1,
'reference_fa' => '/data/public/reference_genomes/genome_hg18.fa',
'dbsnp_rod' => '/data/nextgen1/pipeline/dbsnp_130_hg18.rod',
'target_bed' => '/data/nextgen1/pipeline/targets/SureSelect_All_Exon_G3362_with_names.v2.bed', #Agilent SureSelect Target Enrichment
'picard_dir' => '/home/gmslz/tools/picard-tools-1.44',
'gatk_dir' => '/data/public/tools/gatk-v4333/Sting',
'tmp_dir' => '/home/tmp',
'java_mem' => '-Xmx4g');

$NGS_obj -> check_fastq_format();
#$NGS_obj -> illumina2sanger();
#$NGS_obj -> qcstats();
#$NGS_obj -> fastq_trimmer(1, 74);
#$NGS_obj -> aln_bwa(); #alignment step, which includes bwa, sort, rm_pcrdup, add_rg
$NGS_obj -> snp_analysis_gatk(); #gatk step, which includes realignment, recalibration, genotyping (indel, snp)

NGS Illumina pipeline perl package

#!/usr/bin/perl
use lib "./";
################################################################
# NGS_Illumina Package                                         #
# Author:   Zhengdeng Lei                           #
# E-mail:   leizhengdeng@hotmail.com                           #
################################################################

# Example in caller NGS_Illumina.pl
# Assume the paths for fastx_tk, maq, bwa, samtools are in ENV{PATH}; check $PATH setting $HOME/.bash_profile
# Assume perl and python are installed
#my $NGS_obj = new NGS_Illumina ('sample_name' => 'Testis_T28',
# 'input_fq1' => 'Run1_testicular-28T_lane2_read1_sequence.txt',
# 'input_fq2' => 'Run1_testicular-28T_lane2_read2_sequence.txt',
# 'pe' => 1,
# 'reference_fa' => '/data/public/reference_genomes/genome_hg18.fa',
# 'dbsnp_rod' => '/data/nextgen1/pipeline/dbsnp_130_hg18.rod',
# 'picard_dir' => '/home/gmslz/tools/picard-tools-1.44',
# 'gatk_dir' => '/data/public/tools/gatk-v4333/Sting',
# 'tmp_dir' => '/home/tmp',
# 'java_mem' => '-Xmx4g');
#
#$NGS_obj -> check_fastq_format();
#$NGS_obj -> illumina2sanger();
#$NGS_obj -> qcstats();
#$NGS_obj -> fastq_trimmer(1, 74);
#$NGS_obj -> aln_bwa(); #alignment step, which includes bwa, sort, rm_pcrdup, add_rg
#$NGS_obj -> snp_analysis_gatk(); #gatk step, which includes realignment, recalibration, genotyping (indel, snp)

package NGS_Illumina;
#use strict;
use warnings;
use Time::Local;
use Time::localtime;

# private members:
# sample_name <sample name>
# input_fq1 <input fastq1>
# input_fq2 <input fastq2>
# sanger_fq1 <input fastq1 in sanger format>
# sanger_fq2 <input fastq1 in sanger format>
# reference_fa <reference genome>
# fastq_format <fastq format>
# "Sanger" -- quality score from 0 to 93 using ASCII 33 to 126
# "Solexa" Illumina 1.0 -- quality score from -5 to 62 using ASCII 59 to 126
# "Illumina" Illumina 1.3+ -- quality score from 0 to 62 using ASCII 64 to 126
# pe <paired end>
# 1 paired end
# 0 single end
#my @attributes = qw/reference_fa sample_name input_fq1 input_fq2 fastq_format sanger_fq1 sanger_fq2 pe/;

######################################################################################################
# The subroutines "new" and "initialize" for the constructor
sub new  
{
my $class_name = shift; # Gets class name from parmlist
my $this = {}; # Creates reference to object
bless($this, $class_name); # Associates reference with class
$this->initialize(@_); # Call local initialize subroutine
# passing rest of parmlist
return $this; # Explicit return of value of $this
}

sub initialize
{
my $this = shift;    # Receive an object reference as 1st param
my %parms = @_;         # Get passed parameters from the call to "new"

# Change hash reference for key with supplied value, or assign default
$this->{reference_fa} = $parms{reference_fa} || 'hg18.fa';
$this->{sample_name} = $parms{sample_name} || 'undef_sample_name';
$this->{input_fq1} = $parms{input_fq1} || 0;
$this->{input_fq2} = $parms{input_fq2} || 0;
$this->{pe} = $parms{pe} || 1;
$this->{fastq_format} = $parms{fastq_format} || 'Illumina';
$this->{sanger_fq1} = $parms{sample_name}.'_read1_sanger.fq';
$this->{sanger_fq2} = $parms{sample_name}.'_read2_sanger.fq';
$this->{dbsnp_rod} = $parms{dbsnp_rod} || '/data/nextgen1/pipeline/dbsnp_130_hg18.rod';
$this->{target_bed} = $parms{target_bed} || '/data/nextgen1/pipeline/targets/SureSelect_All_Exon_G3362_with_names.v2.bed';

$this->{java_mem} = $parms{java_mem} || '-Xmx4g';
$this->{num_thread} = $parms{num_thread} || 8;
$this->{gatk_dir} = $parms{gatk_dir} || '/data/public/tools/gatk-v4333/Sting';
$this->{picard_dir} = $parms{picard_dir} || '/home/gmslz/tools/picard-tools-1.44';
$this->{tools_dir} = $parms{tools_dir} || '/data/public/tools/gatk_pipeline/src-pipeline';
$this->{tmp_dir} = $parms{tmp_dir} || '/home/tmp';

$this->{picard_MarkDuplicates} = "java $this->{java_mem} -jar $this->{picard_dir}/MarkDuplicates.jar";
$this->{picard_AddOrReplaceReadGroups} = "java $this->{java_mem} -jar $this->{picard_dir}/AddOrReplaceReadGroups.jar";
$this->{picard_SortSam} = "java $this->{java_mem} -jar $this->{picard_dir}/SortSam.jar";
$this->{picard_FixMateInformation} = "java $this->{java_mem} -Djava.io.tmpdir=$this->{tmp_dir}/ -jar $this->{picard_dir}/FixMateInformation.jar";
$this->{gatk_GenomeAnalysisTK} = "java $this->{java_mem} -Djava.io.tmpdir=$this->{tmp_dir}/ -jar $this->{gatk_dir}/dist/GenomeAnalysisTK.jar";
$this->{gatk_makeIndelMask_py} = "$this->{gatk_dir}/python/makeIndelMask.py";

$this->{aligned_bam_file} = "./bam/$this->{sample_name}.sorted.rmdup.rg.bam"; #This is the final bam file after bwa alignment, i.e. before gatk
$this->{final_recalibrated_bam_file} = "./bam/$this->{sample_name}.final.recalibrated.bam"; #This is the final recalibrated bam file after after gatk

if ($this->{fastq_format} eq "Sanger")
{
$this->{sanger_fq1} = $this->{input_fq1};
$this->{sanger_fq2} = $this->{input_fq2};
}


#print $ENV{'PATH'}, "\n\n\n";

my @subfolders = qw/qcstats bam snp_analysis_gatk/;
foreach my $dir (@subfolders)
{
if (!(-e $dir))
{
mkdir($dir);
}
}
}

sub modify  #  modify the object
{
my $this = shift;  # Figure out who I am
my %parms = @_;       # Make hash out of rest of parm list
for(keys %parms)
{
$this->{$_} = $parms{$_}; # Replace with new value
}
}

sub snp_analysis_gatk
{
# GATK Documents:
# http://www.broadinstitute.org/gatk/gatkdocs/
# e.g.:
# http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_filters_DuplicateReadFilter.html
# http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_filters_VariantFiltration.html

my $this = shift;
#chdir("../snp_analysis_gtk");
my $max_depth = 10000;

# goto EVAL1;
run_cmd("$this->{gatk_GenomeAnalysisTK} -T DepthOfCoverage -l INFO ". #not support -nt $this->{num_thread}
"-D $this->{dbsnp_rod} -R $this->{reference_fa} -L $this->{target_bed} ".
"-I $this->{aligned_bam_file} ".
"-o ./snp_analysis_gatk/$this->{sample_name}.depthofcoverage > ./log/$this->{sample_name}.depthofcoverage.log");







#step 1 --- Local realignment around indels
#step 1.1 --- gatk -T RealignerTargetCreator; Determining (small) suspicious intervals around indels
run_cmd("$this->{gatk_GenomeAnalysisTK} -T RealignerTargetCreator -l INFO ". #not support -nt $this->{num_thread}
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
"-I $this->{aligned_bam_file} ".
"-o ./snp_analysis_gatk/$this->{sample_name}.intervals > ./log/$this->{sample_name}.intervals.log");
#step 1.2 --- gatk -T IndelRealigner; Running the realigner over those intervals
run_cmd("$this->{gatk_GenomeAnalysisTK} -T IndelRealigner -l INFO ". #not support -nt $this->{num_thread}
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
"-I $this->{aligned_bam_file} -targetIntervals ./snp_analysis_gatk/$this->{sample_name}.intervals ".
"-o ./bam/$this->{sample_name}.realigned.bam -stats ./bam/$this->{sample_name}.realign.stats.txt > ./log/$this->{sample_name}.realign.log");
run_cmd("samtools index ./bam/$this->{sample_name}.realigned.bam");
#step 1.3 --- picard SortSam
run_cmd("$this->{picard_SortSam} SO=coordinate VALIDATION_STRINGENCY=SILENT ".
"I=./bam/$this->{sample_name}.realigned.bam ".
"O=./bam/$this->{sample_name}.realigned.sorted.bam");
run_cmd("samtools index ./bam/$this->{sample_name}.realigned.sorted.bam");
#step 1.4 --- picard FixMateInformation; Fixing the mate pairs of realigned reads
run_cmd("$this->{picard_FixMateInformation} SO=coordinate VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=2000000 ".
"I=./bam/$this->{sample_name}.realigned.sorted.bam ".
"O=./bam/$this->{sample_name}.realigned.sorted.fixed.bam");
run_cmd("samtools index ./bam/$this->{sample_name}.realigned.sorted.fixed.bam");
#step 1.5 --- skip picard MarkDuplicates


#step 2 --- Recalibrate quality scores --> Final Recalibrated Bam
#step 2.1 --- gatk -T CountCovariates;
run_cmd("$this->{gatk_GenomeAnalysisTK} -T CountCovariates -l INFO -nt $this->{num_thread} ".
"--downsample_to_coverage $max_depth --default_platform Illumina ".
"-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate ".
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
  "-I ./bam/$this->{sample_name}.realigned.sorted.fixed.bam ".
"-recalFile ./snp_analysis_gatk/$this->{sample_name}.recal.csv > ./log/$this->{sample_name}.count.covariates.log");
#step 2.2 --- gatk -T TableRecalibration  --> Final Recalibrated Bam
run_cmd("$this->{gatk_GenomeAnalysisTK} -T TableRecalibration -l INFO --default_platform Illumina ". #not support -nt $this->{num_thread}
"-R $this->{reference_fa} ".
  "-I ./bam/$this->{sample_name}.realigned.sorted.fixed.bam -recalFile ./snp_analysis_gatk/$this->{sample_name}.recal.csv ".
"--out $this->{final_recalibrated_bam_file} > ./log/$this->{sample_name}.recalibration.log");
run_cmd("samtools index $this->{final_recalibrated_bam_file}");
# skip 2.3 -T CountCovariates
# skip 2.4 AnalyzeCovariate

#############################################################################
# You may merge multiple final_recalibrated_bam files before Genotyping #
#############################################################################


#step 3 --- Genotyping
#step 3.1 --- gatk -T IndelGenotyperV2;
#print STDERR "\nRunning the indel genotyper:\n";
    run_cmd("$this->{gatk_GenomeAnalysisTK} -T IndelGenotyperV2 -l INFO --window_size 450 ". #not suport -nt $this->{num_thread}
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
"-I $this->{final_recalibrated_bam_file} ".
"-o ./snp_analysis_gatk/$this->{sample_name}.indels.vcf -bed ./snp_analysis_gatk/$this->{sample_name}.indels.bed > ./log/$this->{sample_name}.indels.log");

# skip filter indels: perl /tools/filterSingleSampleCalls.pl
#step 3.2 python/makeIndelMask.py
# create intervals to mask snvs close to indels
# The output of the IndelGenotyper is used to mask out SNPs near indels.
# The number 10 in this command stands for the number of bases that will be included on either side of the indel
run_cmd("python $this->{gatk_dir}/python/makeIndelMask.py ./snp_analysis_gatk/$this->{sample_name}.indels.bed 10 ./snp_analysis_gatk/$this->{sample_name}.indels.mask.bed");
#step 3.3 --- gatk -T Unifiedgenotyper;
run_cmd("$this->{gatk_GenomeAnalysisTK} -T UnifiedGenotyper -l INFO ". #not suport -nt $this->{num_thread}
"-stand_call_conf 30 -stand_emit_conf 10 -pl SOLEXA -mmq 30 -mm40 3 ".
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
"-I $this->{final_recalibrated_bam_file} ".
"-o ./snp_analysis_gatk/$this->{sample_name}.raw.snps.vcf > ./log/$this->{sample_name}.snp.log");

#step 3.4 --- gatk -T VariantFiltration; VariantFiltration is used to annotate suspicious calls from VCF files based on their failing given filters.
run_cmd("$this->{gatk_GenomeAnalysisTK} -T VariantFiltration -l INFO ". #not suport -nt $this->{num_thread}
'--clusterSize 3 --clusterWindowSize 10 '.
'--filterExpression "DP <= 8" --filterName "DP8" '.
'--filterExpression "SB > -0.10" --filterName "StrandBias" '.
'--filterExpression "HRun > 8" --filterName "HRun8" '.
'--filterExpression "QD < 5.0" --filterName "QD5" '.
'--filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filterName "hard2validate" '.
"-D $this->{dbsnp_rod} -R $this->{reference_fa} ".
"-B:variant,VCF ./snp_analysis_gatk/$this->{sample_name}.raw.snps.vcf ".
"-B:mask,Bed ./snp_analysis_gatk/$this->{sample_name}.indels.mask.bed --maskName InDel ".
"-o ./snp_analysis_gatk/$this->{sample_name}.snp.filtered.vcf > ./log/$this->{sample_name}.snp.filter.log");

## grep -i 'indel'  Testis_T28.snp.filtered.vcf|more
## grep -i 'pass'  Testis_T28.snp.filtered.vcf|more

#EVAL1:
run_cmd("$this->{gatk_GenomeAnalysisTK} -T VariantEval -l INFO ". #not suport -nt $this->{num_thread}
"-D $this->{dbsnp_rod} -R $this->{reference_fa} -L $this->{target_bed} ".
#-G Indicates that your VCF file has genotypes; -A Print extended evaluation information; -V Print a list of interesting sites (FPs, FNs, etc).
"-B:eval,VCF ./snp_analysis_gatk/$this->{sample_name}.snp.filtered.vcf ".
"-o ./snp_analysis_gatk/$this->{sample_name}.snp.filtered.eval > ./log/$this->{sample_name}.snp.filter.eval.log");



#annotate variants Annovar
#$cmd = "annotate_variation.pl -filter -batchsize 50m -dbtype snp$verdbsnp -buildver $buildver -outfile $outfile $queryfile $dbloc";


# Consensus Quality: Phred-scaled likelihood that the genotype is wrong
# SNP Quality: Phred-scaled likelihood that the genotype is identical to the reference,
# which is also called `SNP quality'. Suppose the reference base is A and in alignment
# we see 17 G and 3 A. We will get a low consensus quality because it is difficult to
# distinguish an A/G heterozygote from a G/G homozygote. We will get a high SNP
# quality, though, because the evidence of a SNP is very strong.
# RMS: root mean square (RMS) mapping quality, a measure of the variance of quality scores
# Coverage: reads covering the position
}


########################################################################################################
# 0. check fastq format
sub check_fastq_format
{
my $this = shift;  # Figure out who I am
open(FASTQ, $this->{input_fq1}) or die "can't open $this->{input_fq1}: $@\n";

my $fastq_format;
    my $counter = 2000; # count 2000 lines of quality scores
my $min_qualscore = 1e6;
my $max_qualscore = -1e6;
    my $tmp;
    while ($counter > 0)
{
$tmp = <FASTQ>; # @read name
$tmp = <FASTQ>; # base calls
$tmp = <FASTQ>; # +read name
my $scores = <FASTQ>; # quality scores
if (!$scores) { last }
#print $scores;
chomp($scores);
for my $chr (map(ord, split(//, $scores)))
{
if ($chr < $min_qualscore)
{
$min_qualscore = $chr;
}
if ($chr > $max_qualscore)
{
$max_qualscore = $chr;
}

}
$counter--;
}
# Phred+33 means quality score + 33 = ASCII
if ($min_qualscore >= 33 && $max_qualscore <= 126)
{
$fastq_format = "Sanger";
}

# Solexa+64 means quality score + 64 = ASCII
if ($min_qualscore >= 59 && $max_qualscore <= 126)
{
$fastq_format = "Solexa";
}

# Phred+64 means quality score + 64 = ASCII
if ($min_qualscore >= 64 && $max_qualscore <= 126)
{
$fastq_format = "Illumina 1.3+";
}

# Phred+64 in both both Illumina 1.3+ and Illumina 1.5+
#if ($min_qualscore >= 66 && $max_qualscore <= 126)
#{
# $fastq_format = "Illumina 1.5+";
#}

# Illumina 1.8+ returned to the use of the Sanger format (Phred+33)
print "$fastq_format fastq format: ASCII(min, max) = ($min_qualscore, $max_qualscore)\n";
$this->{fastq_format} = $fastq_format;
if ($this->{fastq_format} eq "Sanger")
{
$this->{sanger_fq1} = $this->{input_fq1};
$this->{sanger_fq2} = $this->{input_fq2};
}
return $fastq_format;
}

########################################################################################################
# 1. QCGroomer step: illumina2sanger
# Alternative converter "fastq_quality_converter -i $this->{input_fq1} -o $this->{sanger_fq1}"
sub illumina2sanger
{
my $this = shift;  # Figure out who I am
my $fastq_converter;
if ($this->{fastq_format} eq "Sanger")
{
print "The input fastq is already in Sanger format, no need for conversion!!\n";
return;
}
elsif ($this->{fastq_format} eq "Solexa")
{
$fastq_converter = "maq sol2sanger";
}
elsif ($this->{fastq_format} eq "Illumina 1.3+")
{
$fastq_converter = "maq ill2sanger";
}
else
{
print "No Sanger or Illumina fastq format found, exit now!!\n";
exit(0);
}
run_cmd("$fastq_converter $this->{input_fq1} $this->{sanger_fq1}");

if($this->{pe})
{
run_cmd("$fastq_converter $this->{input_fq2} $this->{sanger_fq2}");
}

}


########################################################################################################
# 2. QCStats step: qcstats
# use fastx tool kit to process
# http://hannonlab.cshl.edu/fastx_toolkit/commandline.html
# check linux version by "uname -a"
# To set PATH=$PATH:$HOME/bin:$HOME/tools/fastx_tk
# vi .bash_profile
# source .bash_profile
sub qcstats
{
my $this = shift;
run_cmd("fastx_quality_stats -Q 33 -i $this->{sanger_fq1} -o ./qcstats/$this->{sanger_fq1}"."_stats.txt"); #option -Q 33 is for Sanger format
run_cmd("fastq_quality_boxplot_graph.sh -i ./qcstats/$this->{sanger_fq1}_stats.txt -o ./qcstats/$this->{sanger_fq1}_stats.png -t \"$this->{sanger_fq1}\"");
# my $qcnecleotide_cmd1 = "fastx_nucleotide_distribution_graph.sh -i $this->{sanger_fq1}"."_stats.txt -o ./qcstats/$this->{sanger_fq1}_nuc.png -t \"$this->{sanger_fq1} nucleotide distribution\"";
if($this->{pe})
{
run_cmd("fastx_quality_stats -Q 33 -i $this->{sanger_fq2} -o ./qcstats/$this->{sanger_fq2}"."_stats.txt"); #option -Q 33 is for Sanger format
run_cmd("fastq_quality_boxplot_graph.sh -i ./qcstats/$this->{sanger_fq2}_stats.txt -o ./qcstats/$this->{sanger_fq2}_stats.png -t \"$this->{sanger_fq2}\"");
# my $qcnecleotide_cmd2 = "fastx_nucleotide_distribution_graph.sh -i $this->{sanger_fq2}"."_stats.txt -o ./qcstats/$this->{sanger_fq2}_nuc.png -t \"$this->{sanger_fq2} nucleotide distribution\"";
}
}

sub fastq_trimmer
{
my $this = shift;
my ($first_base_pos, $last_base_pos) = @_;
print "reads at positions($first_base_pos, $last_base_pos) selected\n";
my $trimmed_fq1 = $this->{sanger_fq1};
$trimmed_fq1 =~ s/fq$/trimmed.fq/;
run_cmd("fastx_trimmer -Q 33 -f $first_base_pos -l $last_base_pos -i $this->{sanger_fq1} -o $trimmed_fq1");
$this->{sanger_fq1} = $trimmed_fq1;

if($this->{pe})
{
my $trimmed_fq2 = $this->{sanger_fq2};
$trimmed_fq2 =~ s/fq$/trimmed.fq/;
run_cmd("fastx_trimmer -Q 33 -f $first_base_pos -l $last_base_pos -i $this->{sanger_fq2} -o $trimmed_fq2");
$this->{sanger_fq2} = $trimmed_fq2;
}
}

sub aln_bwa
{
#    my (%sample_setting) = @_;
my $this = shift;
run_cmd("bwa aln -t $this->{num_thread} $this->{reference_fa} $this->{sanger_fq1} > ./bam/$this->{sample_name}.read1.sai");

if($this->{pe})
{
run_cmd("bwa aln -t $this->{num_thread} $this->{reference_fa} $this->{sanger_fq2} > ./bam/$this->{sample_name}.read2.sai");
#to add read group information -r '@RG\tID:foo\tSM:bar'
run_cmd("bwa sampe $this->{reference_fa} ./bam/$this->{sample_name}.read1.sai ./bam/$this->{sample_name}.read2.sai ".
"$this->{sanger_fq1} $this->{sanger_fq2} > ./bam/$this->{sample_name}.sam");
}
else
{
run_cmd("bwa samse $this->{reference_fa} ./bam/$this->{sample_name}.read1.sai $this->{sanger_fq1} > ./bam/$this->{sample_name}.sam");
}
# #//sed 's/[ \t]*$//g' $sample_name.sam >$sample_name.sed.sam
run_cmd("samtools import $this->{reference_fa}.fai ./bam/$this->{sample_name}.sam ./bam/$this->{sample_name}.bam");
run_cmd("samtools sort ./bam/$this->{sample_name}.bam ./bam/$this->{sample_name}.sorted");
run_cmd("samtools index ./bam/$this->{sample_name}.sorted.bam");

#scp -r ./picard-tools-1.44/ NUSSTF\\gmslz@172.25.138.143:~/tools/
#add $HOME/tools/picard-tools-1.44 to PATH in file .bash_profile
#source .bash_profile
##### remove PCR duplicate using PICARD
##### remove PRC duplicates by picard, use AS=true so that picard knows the bam is sorted.
run_cmd("$this->{picard_MarkDuplicates} VALIDATION_STRINGENCY=SILENT M=PCR_duplicates.M.txt REMOVE_DUPLICATES=true AS=true ".
"I=./bam/$this->{sample_name}.sorted.bam ".
"O=./bam/$this->{sample_name}.sorted.rmdup.bam");

#add read group inforation
run_cmd("$this->{picard_AddOrReplaceReadGroups} RGID=lane1 RGLB=$this->{sample_name}.fastq RGPL=Illumina RGPU=run RGSM=$this->{sample_name} CREATE_INDEX=TRUE VALIDATION_STRINGENCY=SILENT ".
"I=./bam/$this->{sample_name}.sorted.rmdup.bam ".
"O=./bam/$this->{sample_name}.sorted.rmdup.rg.bam");
run_cmd("samtools index ./bam/$this->{sample_name}.sorted.rmdup.rg.bam");

unlink("./bam/$this->{sample_name}.read1.sai", "./bam/$this->{sample_name}.read2.sai", "./bam/$this->{sample_name}.sam");
}



sub run_cmd
{
my $cmd = shift;
my $start_time = localtime();
my ($year, $mon, $day, $hour, $min, $sec) = ($start_time->year + 1900, $start_time->mon+1, $start_time->mday, $start_time->hour, $start_time->min, $start_time->sec);
my $start_time_sec = timelocal($sec, $min, $hour, $day, $mon, $year);
printf("%04d/%02d/%02d %02d:%02d:%02d", $year, $mon, $day, $hour, $min, $sec);
print "\tSTART $cmd\n";

my $cmd_res = system($cmd);
if ($cmd_res)
{
my $error_time = localtime;
($year, $mon, $day, $hour, $min, $sec) = ($error_time->year + 1900, $error_time->mon+1, $error_time->mday, $error_time->hour, $error_time->min, $error_time->sec);
printf("%04d/%02d/%02d %02d:%02d:%02d", $year, $mon, $day, $hour, $min, $sec);
print "\tERROR: $?\nEXIT!!\n\n";
exit 1;

}
else
{
my $end_time = localtime;
($year, $mon, $day, $hour, $min, $sec) = ($end_time->year + 1900, $end_time->mon+1, $end_time->mday, $end_time->hour, $end_time->min, $end_time->sec);
my $end_time_sec = timelocal($sec, $min, $hour, $day, $mon, $year);
my $difference = $end_time_sec - $start_time_sec;
my $seconds    =  $difference % 60;
$difference = ($difference - $seconds) / 60;
my $minutes    =  $difference % 60;
$difference = ($difference - $minutes) / 60;
my $hours      =  $difference % 24;
printf("%04d/%02d/%02d %02d:%02d:%02d", $year, $mon, $day, $hour, $min, $sec);
print "\tSUCCESS after running $hours hours $minutes minutes $seconds seconds\n\n";
}
}

1;


#
# Reference:
# http://www.bbmriwiki.nl/wiki/PipelineCommands
# http://www.bbmriwiki.nl/wiki/SnpCallingPipeline/VariantCalling
# http://www.broadinstitute.org/gatk/gatkdocs/

Mammary development meets cancer genomics

Mammary development meets cancer genomics

Build NTP predictor using order statistics

To output of limma
1. (for universal platform) convert to ENTREZ GENE ID, remove those with opposite direction, i.e. opposite t-score for same gene
2. calculate ranking.p by using p.val and abs(FC) based on order statistics
3. using -log10(ranking.p) as feature weight.