外显子测序数据分析流程

1、数据下载
2、与参考序列比对（采用bwa）
3、标记PCR重复
4、碱基矫正
5、质量控制
6、VEP标记

数据下载

数据下载方式：

/mnt/raid61/Personal_data/zhouran/soft/linuxnd login -u X101SC20030543-Z01-J077 -p cxpez225 /mnt/raid61/Personal_data/zhouran/soft/linuxnd cp -d oss://CP2018091214925/H101SC20030543/RSCS7000/X101SC20030543-Z01/X101SC20030543-Z01-J077/ /mnt/raid62/BetaCoV/Person/zengwanqin/datalss

1、比对

ls /mnt/raid62/BetaCoV/Person/zengwanqin/datalss/X101SC20030543-Z01-J077/2.cleandata/夏天翼2号管_FDHE202653854-1a/ 夏天翼2号管_FDHE202653854-1a_1.clean.fq.gz>1 ls /mnt/raid62/BetaCoV/Person/zengwanqin/datalss/X101SC20030543-Z01-J077/2.cleandata/夏天翼2号管_FDHE202653854-1a/ 夏天翼2号管_FDHE202653854-1a_2.clean.fq.gz>2 cut -d"/" -f 7 1 |cut -d"_" -f 1 > 0 paste 0 1 2 > config source activate wes INDEX=/mnt/raid62/Personal_data/zhangyiming/GATK/hg38/Homo_sapiens_assembly38.fasta cat config |while read id do arr=($id) fq1${arr[1]} fq2=${arr[2]} sample=${arr[0]} echo $sample $fq1 $fq2 bwa mem -t 5 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2 | samtools sort -@ 5 -o $sample.bam - done

##质量控制

time $gatk CollectHsMetrics \ -BI /mnt/raid64/pancreas_sc/align/WES/info/Agilent/S07604514_Covered.interval_list \ -I FDHE202653854.bam \ -O xiatianyi.txt.metrics \ -TI /mnt/raid64/pancreas_sc/align/WES/info/Agilent/S07604514_Covered.interval_list \ --VALIDATION_STRINGENCY=LENIENT查看文件：
samtools view 夏天翼2号管.bam | less
samtools view -H 夏天翼2号管.bam |grep -v "SQ"
samtools index 夏天翼2号管.bam

2、标记重复序列####定义变量

ref=/mnt/raid62/Personal_data/zhangyiming/GATK/hg38/Homo_sapiens_assembly38.fasta
snp=/mnt/raid62/Personal_data/zhangyiming/GATK/hg38/dbsnp_146.hg38.vcf.gz
indel=/mnt/raid62/Personal_data/zhangyiming/GATK/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
gatk=/mnt/raid61/Personal_data/zhouran/soft/gatk-4.1.0.0/gatk

标记PCR重复reads

sample=xiatianyi echo $sample time $gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \ -I FDHE202653854.bam \ -O ${sample}_marked.bam \ -M $sample.metrics \ 1>log.mark 2>&1

#FixMateInformation

time $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

#进行index

samtools index ${sample}_marked_fixed.bam

#BaseRecalibrator 碱基矫正

time $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 time $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

#HaplotypeCaller命令

$gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \ -R $ref \ -I ${sample}_bqsr.bam \ --dbsnp $snp \ -O ${sample}_raw.vcf \ 1>${sample}_log.HC 2>&1

#gtf文件：/mnt/raid64/ref_genomes/HomSap/release101/Homo_sapiens.GRCh38.101.sorted.gtf.gz
#VEP进行功能注释

/home/zhouran/data/miniconda3/bin/perl /home/zhouran/vep/vep --species homo_sapiens --assembly GRCh38 --no_progress --no_stats --buffer_size 5000 --sift b --ccds --uniprot --hgvs --symbol --numbers --domains --gene_phenotype --canonical --protein --biotype --uniprot --tsl --variant_class --shift_hgvs 1 --check_existing --total_length --allele_number --no_escape --xref_refseq --failed 1 --vcf --flag_pick_allele --pick_order canonical,tsl,biotype,rank,ccds,length --dir /home/zhouran/.vep --fasta /home/zhouran/.vep/homo_sapiens/100_GRCh38/Homo_sapiens.GRCh38.100.dna.primary_assembly.fa.gz --format vcf --input_file xiatianyi_raw.vcf --output_file xiatianyi_raw.vep.vcf.gz --offline --pubmed --fork 4 --cache_version 100 --polyphen b --af --af_1kg --af_esp --af_gnomad --regulatory

##保留probe内的SNP index

bcftools view xiatianyi_raw.vep.vcf.gz -R /mnt/raid64/pancreas_sc/align/WES/info/Agilent/S07604514_Covered.bed | bgzip > subset.vep.vcf.gz

##注释

vep-annotation-reporter subset.vep.vcf.gz Consequence SYMBOL Feature -o xiatianyi.subset.tsv

##bam-readcount

/mnt/raid62/BetaCoV/Person/zengwanqin/datalss/X101SC20030543-Z01-J077/2.cleandata/FDHE202653854-1a/bam-readcount/bin/bam-readcount \ -f /mnt/raid62/Personal_data/zhangyiming/GATK/hg38/Homo_sapiens_assembly38.fasta \ -l /mnt/raid64/pancreas_sc/align/WES/info/Agilent/S07604514_Covered.bed xiatianyi_bqsr.bam > subset_readcounts.tsv bam-readcount -f /mnt/raid62/Personal_data/zhangyiming/GATK/hg38/Homo_sapiens_assembly38.fasta xiatianyi_bqsr.bam chr1:11000-12000 > read_count.tsv

##将VCF和bam-readcounts合并##https://vatools.readthedocs.io/en/latest/vcf_readcount_annotator.html

vcf-readcount-annotator xiatianyi_raw.vep.vcf subset_readcounts.tsv DNA \ -t all -o all_annotated_vcf bcftools query -f '[%AD\t%DP]\n' all_annotated_vcf > all_annotated_vcf_ad_pd.tsv vep-annotation-reporter all_annotated_vcf Consequence SYMBOL Feature -o xiatianyi.subsetrc.tsv

##用R合并xiatianyi.subsetrc.tsv all_annotated_vcf_ad_pd.tsv提出aa和密码子

vep-annotation-reporter subset.vep.vcf.gz Consequence SYMBOL Feature HGVSp Codons -o xiatianyi.subset.tsv

#看该区域有没有突变
/

home/zhouran/data/anaconda3/bin/tabix xiatianyi_raw.vep.vcf.gz chr1:160288594-160343273| less
ls -lh xiatianyi*.bam|cut -d " " -f 5-samtools view -h FDHE202653854.bam chr1:160288594-160343273|samtools sort -o xiatianyi.copa.bam -
samtools view -h xiatianyi_bqsr.bam chr1:160288594-160343273|samtools sort -o xiatianyi_bqsr.copa.bam -
samtools view -h xiatianyi_marked.bam chr1:160288594-160343273|samtools sort -o xiatianyi7_marked.copa.bam -
samtools view -h xiatianyi_marked_fixed.bam chr1:160288594-160343273|samtools sort -o xiatianyi_marked_fixed.copa.bam -
ls -lh *copa.bam
ls *.copa.bam|xargs -i samtools index {}

zcat /mnt/raid64/ref_genomes/HomSap/release101/Homo_sapiens.GRCh38.101.sorted.gtf.gz |grep -w STAT3|head|less -SN

/home/zhouran/data/anaconda3/bin/tabix xiatianyi_raw.vep.vcf.gz chr17:42313324-42388482| less
zcat /mnt/raid64/ref_genomes/HomSap/release101/Homo_sapiens.GRCh38.101.sorted.gtf.gz |grep -w LRBA|head|less -SN