WES(全外显子)分析(上)

感悟：

软件不看帮助文档，不阅读说明书，就只能抄代码，却不知道错在哪里。不想要成为代码搬运工。

1. 分析前工作准备

1.1 创建环境：

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

#举个栗子
conda search softwarename
conda install -y softwarename

软件安装

2. 数据下载

找SRA数据,下载SRR list

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

nohup bash download. sh &

3. sra转为fastq

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

ls SRR5660416.sra |fastq-dump --gzip --split-3 -O ./ SRR5660416.sra #单个测试
cat >sra.sh #写脚本
ls * |while read id; do (fastq-dump --gzip --split-3 -O ./ ${id} 1>${id}.log 2>&1);done #多个循环
nohup bash sra.sh & #挂后台运行

4. 质控

4.1 fastqc 质量报告

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得到文件(举个栗子)

4.2 去接头

mkdir clean && cd clean
#单个测试
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
#多个循环 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 &

质控，去接头

下述流程，均用一个sample少量序列做分析，故给sample赋值

5. 有参Mapping

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

#从gz到sam到sort到bam，注意输出文件后的-必须加！！！！思考：赋值时的引号，加不加有什么区别？？
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统计
#单个测试
samtools flagstat /home/qmcui/7E5241.chr1_2.sort.bam
#多个
ls *.bam | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done
结果中mapped

6. mark pcr重复

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 标记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

最后编辑于：2019.01.25 15:57:49