感悟：

軟件不看幫助文檔，不閱讀說明書，就只能抄代碼，卻不知道錯(cuò)在哪里。不想要成為代碼搬運(yùn)工。

1. 分析前工作準(zhǔn)備

1.1 創(chuàng)建環(huán)境：

conda create -n wes python=2
conda info --envs
source activate wes

1.2 需要軟件：

ascp,samtools,vcftools,bcftools,fastqc,multiqc,trim_galore,bwa,sra-tools,cutadapt,gatk4

#舉個(gè)栗子
conda search softwarename
conda install -y softwarename

軟件安裝

2. 數(shù)據(jù)下載

找SRA數(shù)據(jù),下載SRR list

cat >download.sh
cat SRR_Acc_List.txt|while read id; do (prefetch ${id}) ;done

nohup bash download. sh &

3. sra轉(zhuǎn)為fastq

mkdir sra && cd sra
mv /home/vip11/ncbi/project/wes/sra/SRR* ./

ls SRR5660416.sra |fastq-dump --gzip --split-3 -O ./ SRR5660416.sra #單個(gè)測試
cat >sra.sh #寫腳本
ls * |while read id; do (fastq-dump --gzip --split-3 -O ./ ${id} 1>${id}.log 2>&1);done #多個(gè)循環(huán)
nohup bash sra.sh & #掛后臺(tái)運(yùn)行

4. 質(zhì)控

4.1 fastqc 質(zhì)量報(bào)告

mkdir tmp && cd tmp
fastqc /home/qmcui/7E5240_L1_A001.L1_1.fastq.gz /home/qmcui/7E5240_L1_A001.L1_2.fastq.gz -o ./tmp

fastqc得到文件(舉個(gè)栗子)

4.2 去接頭

mkdir clean && cd clean
#單個(gè)測試
trim_galore --phred33 -q 25 -e 0.1 --length 36 --stringency 3 --paired -o ~/ ../sra/tmp/7E5241.L1_1.fastq.gz ../sra/tmp/7E5241.L1_2.fastq.gz
#多個(gè)循環(huán) trim_galore
ls /home/qmcui/*1.fastq.gz>./1;cat 1
ls /home/qmcui/*2.fastq.gz>./2;cat 2
paste 1 2 >config
cat >trim.sh
cat config|while read id;do arr=(${id}); fq1=${arr[0]}; fq2=${arr[1]}; echo $fq1 $fq2 trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o ../clean  $fq1 $fq2 >./${id}.log 2>&1;done
nohup bash mvadaptor.sh &

質(zhì)控，去接頭

下述流程，均用一個(gè)sample少量序列做分析，故給sample賦值

5. 有參Mapping

mkdir sm_sort_bam && cd sm_sort_bam
#先建立好索引，怎么做？？？

#從gz到sam到sort到bam，注意輸出文件后的-必須加?。。。∷伎迹嘿x值時(shí)的引號(hào)，加不加有什么區(qū)別？？
sample=“7E5239”
INDEX=“/home/qmcui/database/reference/index/bwa/hg38”
fq1=“/home/qmcui/7E5239.L1_1_val_1.fq.gz”
fq2=“/home/qmcui/7E5239.L1_2_val_2.fq.gz”
bwa mem -t 1 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2  | samtools sort -@ 1 -o $sample.bam -

bam統(tǒng)計(jì)
#單個(gè)測試
samtools flagstat /home/qmcui/7E5241.chr1_2.sort.bam
#多個(gè)
ls *.bam | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done
結(jié)果中mapped

6. mark pcr重復(fù)

mkdir mark_bam && cd mark_bam
sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates -I /home/qmcui/7E5241.chr1_2.sort.bam -O ${sample}_marked.bam -M ${sample}.metrics 1>${sample}_log.mark 2>&1

7. GATK4 找變異

7.1 標(biāo)記flagstat

mkdir gatk_bam && cd gatk_bam
sample="7E5241"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation  -I /home/qmcui/7E5241.chr1_2.mk.bam  -O  ${sample}_marked_fixed.bam  -SO coordinate 1>${sample}_log.fix 2>&1

7.2 找變異

sample="7E5241" 
ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"
indel="/home/qmcui/tmp/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
snp="/home/qmcui/dbsnp_146.hg38.vcf.gz"
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator -R $ref -I 7E5241_marked_fixed.bam --known-sites $snp --known-sites $indel -O ${sample}_recal.table

7.3 矯正bam

sample="7E5241" 
ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"
nohup gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR -R $ref -I 7E5241_marked_fixed.bam -bqsr 7E5241_recal.table -O ${sample}_bqsr.bam 1>${sample}_log.ApplyBQSR 2>&1 &

7.4

calling variation 得到vcf

mkdir gatk_single_vcf && cd gatk_single_vcf
ref="/home/qmcui/database/reference/index/hisat/hg38/hg38.fa"
"snp="/home/qmcui/dbsnp_146.hg38.vcf.gz
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller -R $ref -I 7E5241_bqsr.bam --dbsnp $snp -O ${sample}_raw.vcf \ 1>${sample}_log.HC 2>&1

Screen Shot 2019-01-21 at 21.36.02.png