GENOME='/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta'
bam='ESCC13-T1_recal.bam'
samtools mpileup -f $GENOME $bam |\
perl -alne '{$pos=int($F[1]/200000); $key="$F[0]\t$pos";$GC{$key}++ if $F[2]=~/[GC]/;$counts_sum{$key}+=$F[3];$number{$key}++;}END{print "$_\t$number{$_}\t$GC{$_}\t$counts_sum{$_}" foreach keys %number}' |\
sort -k1,1 -k 2,2n >T1.windows

normal='/home/jianmingzeng/data/public/escc/bam/ESCC13-N_recal.bam'
GENOME='/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta'
for tumor in /home/jianmingzeng/data/public/escc/bam/*-T*.bam
do
start=$(date +%s.%N)
echo  `date`

file=$(basename $tumor)
sample=${file%%.*} 
echo calling $sample 

samtools mpileup -q 1 -f $GENOME $normal $tumor |\
awk -F"\t" '$4 > 0 && $7 > 0' |\
java -jar ~/biosoft/VarScan/VarScan.v2.3.9.jar copynumber - $sample --mpileup 1

echo  `date`
dur=$(echo "$(date +%s.%N) - $start" | bc)
printf "Execution time for $sample : %.6f seconds" $dur
echo
done 

java -jar ~/biosoft/VarScan/VarScan.v2.3.9.jar  copyCaller  ESCC13-T1_recal.copynumber --output-file varScan.T1.cnv.called

cd ~/biosoft
# http://conifer.sourceforge.net/download.html
mkdir conifer &&  cd conifer 
wget https://downloads.sourceforge.net/project/conifer/CoNIFER%200.2.2/conifer_v0.2.2.tar.gz 
tar zxvf conifer_v0.2.2.tar.gz  
mkdir test && cd test 
wget http://sourceforge.net/projects/conifer/files/sampledata.tar.gz
tar zxvf sampledata.tar.gz
wget http://sourceforge.net/projects/conifer/files/probes.txt 

python -c "import numpy"
python -c "import tables"
# pip install pysam --user
python -c "import pysam"
python -c "import matplotlib"

bed=exon_probe.hg38.gene.bed
for bam in  *bam
do
file=$(basename $bam )
sample=${file%%.*} 
echo $sample
#export total_reads=$(samtools view -F 0x4 $bam | cut -f 1 | sort | uniq | wc -l) 
export total_reads=$(samtools idxstats  $bam|awk -F '\t' '{s+=$3}END{print s}') 
echo The number of reads is $total_reads 
bedtools multicov -bams  $bam -bed $bed |perl -alne '{$len=$F[2]-$F[1];if($len <1 ){print "$.\t$F[4]\t0" }else{$rpkm=(1000000000*$F[4]/($len* $ENV{total_reads}));print "$.\t$F[4]\t$rpkm"}}' >$sample.rpkm.txt
done 

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.20160908.txt
cat CCDS.20160908.txt |perl -alne '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort -i > exon_probe.grch38.gene.bed
sed 's/^/chr/g' exon_probe.grch38.gene.bed >exon_probe.hg38.gene.bed

 cd ~/biosoft/conifer/test/
 python ~/biosoft/conifer/conifer_v0.2.2/conifer.py analyze \
  	--probes ~/biosoft/conifer/test/probes.txt \
  	--rpkm_dir ~/biosoft/conifer/test/sampledata/RPKM_data/  \
  	--output analysis.hdf5 \
  	--svd 6 \
  	--write_svals singular_values.txt \
  	--plot_scree screeplot.png \
  	--write_sd sd_values.txt

 cd ~/biosoft/conifer/test/
 python ~/biosoft/conifer/conifer_v0.2.2/conifer.py call \
  	--input analysis.hdf5 \
  	--output calls.txt

$ head calls.txt
sampleID	chromosome	start	stop	state
Trio2.fa	chr1	196714946	196759357	del
NA18555	chr1	161483684	161640062	dup
NA15510	chr1	155226464	155261728	dup
Trio3.mo	chr1	149398798	149760173	del
NA19240	chr1	110224430	110254970	dup
Trio3.fa	chr10	46998880	47894049	dup
Trio3.fa	chr10	47948713	48382260	dup
Trio4.fa	chr10	124362280	124370763	dup
NA18956	chr10	124342115	124352111	del

 cd ~/biosoft/conifer/test/
 mkdir export_svdzrpkm/
 python ~/biosoft/conifer/conifer_v0.2.2/conifer.py export \
  	--input analysis.hdf5 \
  	--output ./export_svdzrpkm/
## 需要自定义区域进行可视化
cd ~/biosoft/conifer/test/ 
python ~/biosoft/conifer/conifer_v0.2.2/conifer.py   plot \
  	--input analysis.hdf5 \
  	--region chr#:start-stop \
  	--output image.png \
  	--sample sampleID
 ## 也可以一次性画出所有的图
cd ~/biosoft/conifer/test/
 mkdir call_imgs/
 python ~/biosoft/conifer/conifer_v0.2.2/conifer.py  plotcalls \
  	--input analysis.hdf5 
  	--calls calls.txt 
  	--outputdir ./call_imgs/