#压缩fq/fastq
time gzip -c test.fq > test.fq.gz #单线程
time pigz -k *.fq > *.fq.gz #多线程

#md5sum
md5sum *.gz > md5sum
cat *.md5 > md6sum
diff -y  md5sum md6sum

#判断是Phred33还是Phred64，用awk命令
zcat SRR3081198_trimmed.fq.gz | head -100 | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}'


#fastqc
fastqc -t 4 -o /Volumes/WES/WES_OUTPUT/fastQC *.fastq

#multiqc
multiqc -n wes -o /Volumes/WES/WES_mQC /Volumes/WES/WES_QC2/*_fastqc.zip


#trim_galore
trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o /Users/ning/wes_clean   /Volumes/WES/WES_DATA/WGC101045U_combined_R1.fastq  /Volumes/WES/WES_DATA/WGC101045U_combined_R2.fastq #单线程

ls|grep _R1.fastq.gz > fq1
ls|grep _R2.fastq.gz > fq2
paste fq1 fq2 > config
cat config|head
bash trim.sh #多线程


#fastqc_after_trim，同上
#multiqc_after_trim，同上

#bwa

## 先提取小的fq
find  /Volumes/WES/3_WES_CLEAN  -name "*gz" |grep -v '\._' > fq.txt
cat fq.txt |while read id ;do (cat $id|head -10000 > $(basename $id ".gz"));done

## 然后一个小fq文件比对
sample='WGC111640U'
bwa mem -t 5 -R "@RG\tID:$sample\tSM:$sample\tLB:$sample\tPL:Illumina"  /Users/ning/References/GATK_bundle/hg38/hg38/bwa_index/Homo_sapiens_assembly38.fasta.64 WGC111640U_combined_R1_val_1.fq WGC111640U_combined_R2_val_2.fq  | samtools sort -@ 5 -o WGC111640U.bam -
samtools view -H WGC111640U.bam |grep -v "SQ"

## 或者循环批量比对
ls *1.fq.gz > 1
ls *2.fq.gz > 2
paste 1 2 > config
cat 1|cut -d"_" -f1>0
paste 0 1 2 > config
bash bwa.sh #多线程


#GATK4

ref=/Users/ning/References/GATK_bundle/hg38/hg38/Homo_sapiens_assembly38.fasta
snp=/Users/ning/References/GATK_bundle/hg38/hg38/dbsnp_146.hg38.vcf.gz
indel=/Users/ning/References/GATK_bundle/hg38/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
sample=WGC111640U
echo $sample

##Mark Duplicate Reads
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \
-I $sample.bam -O ${sample}_marked.bam \
-M $sample.metrics \
1>log.mark 2>&1

##Fix Mate Information
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation \
    -I ${sample}_marked.bam \
    -O ${sample}_marked_fixed.bam \
    -SO coordinate \
    1>${sample}_log.fix 2>&1 

samtools index WGC111640U_marked_fixed.bam #进行index

##BaseRecalibrator
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./"  BaseRecalibrator \
    -R $ref  \
    -I ${sample}_marked_fixed.bam  \
    --known-sites $snp \
    --known-sites $indel \
    -O ${sample}_recal.table \
    1>${sample}_log.recal 2>&1
 
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./"   ApplyBQSR \
    -R $ref  \
    -I ${sample}_marked_fixed.bam  \
    -bqsr ${sample}_recal.table \
    -O ${sample}_bqsr.bam \
    1>${sample}_log.ApplyBQSR  2>&1 

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./"  BaseRecalibrator \
    -R $ref  \
    -I ${sample}_bqsr.bam  \
    --known-sites $snp \
    --known-sites $indel \
    -O ${sample}_post_recal.table \
    1>${sample}_log.recal.post 2>&1

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" AnalyzeCovariates \
     -before ${sample}_recal.table \
     -after ${sample}_post_recal.table \
     -plots recalibration_${sample}_plots.pdf

##批量进行
bash gatk4.sh