单样品WES流程bash笔记

需要提前准备的软件：bowtie2 samtools annovar

#! /usr/bin/bash

#单样品WES;

#运行方法:首先mkdir 01raw_data;把原始数据放到该目录下,数据的后缀分别修改成相对应的：_1.fq.gz和_2.fq.gz形式；然后cd ../;nohup sh single_wes.sh &;

mkdir 02clean_data 03align_result 04snv_result 05anno_resutl fastqc_report bowtie2 -p 30 -x /data1/gaojw/refseq/bowtie2_index_hg19/hg19 -1 ./01raw_data/*_1.fq.gz -2 ./01raw_data/*_2.fq.gz -S ./03align_result/out.sam 2>./03align_result/alignment.log cd ./03align_result/ samtools view -bS out.sam >out.bam 2>sam2bam.log

# -b是设置输出BAM格式，-S是设置输入SAM格式

 samtools sort -@ 30 -o out.sorted.bam out.bam 2>bamsort.log

#-@ 30 设置排序和压缩时的线程数量，-o FILE 设置最终排序后的输出文件名

samtools rmdup out.sorted.bam out.rmdup.bam

 #去除其他read

samtools index out.rmdup.bam cd ../04snv_result samtools mpileup -f /data1/gaojw/refseq/hg19.fa -A -d 30000 ../03align_result/out.rmdup.bam | /home/gaojw/miniconda2/bin/varscan mpileup2cns --output-vcf 1 --variants >out.rmdup.varscan.vcf

#-f 输入有索引的fasta参考序列 -A 不跳过异常配对 -d 最大测序深度，过滤掉超深度测序的位点 --output-vcf 1表示输出为vcf格式 --variants 仅报告突变位置（仅mpileup2cns可用）

cd ../05anno_resutl convert2annovar.pl -format vcf4 ../04snv_result/out.rmdup.varscan.vcf >out.avinput

#首先把snp-calling步骤的VCF文件转为annovar软件要求的格式

awk '{if ($1!=chr_no || $2!=Start || $3!=End || $4!=Ref || $5!=Alt)print;chr_no=$1 ; Start=$2 ; End=$3 ; Ref=$4 ; Alt=$5;}' out.avinput > uniq.avinput
table_annovar.pl uniq.avinput /home/lexb/Software_zzx/annovar/humandb -buildver hg19 -out annoresult -remove -protocol refGene,avsnp150,1000g2015aug_all,1000g2015aug_eas,1000g2015aug_sas,1000g2015aug_amr,1000g2015aug_afr,1000g2015aug_eur,gnomad_exome_20190125,dbnsfp30a,intervar_20170202,clinvar_20181225 -operation g,f,f,f,f,f,f,f,f,f,f,f 2>annovar.log &

#注释程序，可一次性完成三种类型的注释

# -buildver hg19 表示使用hg19版本

# -out myanno 表示输出文件的前缀为myanno

# -remove 表示删除注释过程中的临时文件

 # -protocol 表示注释使用的数据库，用逗号隔开，且要注意顺序

# -operation 表示对应顺序的数据库的类型（g代表gene-based、r代表region-based、f代表filter-based），用逗号隔开，注意顺序

wait cp /home/gaojw/WES/config.ini ./ python /home/gaojw/src/InterVar/InterVar-  master/Intervar.py -c ./config.ini 2>Intervar.log awk '{if ($1!=chr_no || $2!=Start || $3!=End || $4!=Ref || $5!=Alt)print;chr_no=$1 ; Start=$2 ; End=$3 ; Ref=$4 ; Alt=$5;}' intervar_result.hg19_multianno.txt.intervar > intervar_result_uniq cut -f 8 uniq.avinput >pos_depth sed -i '1i\depth' pos_depth paste annoresult.hg19_multianno.txt intervar_result_uniq | cut -f 1-13,18,22,27-30,35,36,44,45,55,61,90-94,107-109,124,128 >merge_result paste merge_result pos_depth >merge_result2 bash /home/gaojw/WES/filter_annoresult.sh if [ ! -s ./intervar_result_uniq ];then cd .. bash /home/gaojw/WES/fix_bug.sh fi

#如果intervar_result_uniq为空执行fix bug