从原始测序数据到变异检测的完整分析

代码介绍：
准备双端测序数据
prefetch SRR1770413
fastq-dump --split-files SRR1770413
准备参考基因组
curl -O ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
bgzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa
建立索引
bwa index E.coli_K12_MG1655.fa
samtools faidx E.coli_K12_MG1655.fa
序列比对
bwa mem -M -t 1 -R '@RG\tID:Ecoli\tSM:Ecoli\tLB:Ecoli\tPL:ILLUMINA\tPU:run' E.coli_K12_MG1655.fa SRR1770413_1.fastq SRR1770413_2.fastq > SRR1770413.sam
sam转成bam，并进一步过滤
java -Xmx8G -jar picard.jar SortSam I=SRR1770413.sam O=SRR1770413.bam SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true TMP_DIR=./tmp.363296
java -Xmx8G -jar picard.jar MarkDuplicates I=SRR1770413.bam O=SRR1770413.dup.bam M=SRR1770413.dup.txt REMOVE_DUPLICATES=true VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true TMP_DIR=./tmp.690131
bam文件质量检查
samtools flagstat SRR1770413.dup.bam > SRR1770413.dup.bam.stat
建立dict
java -jar picard.jar CreateSequenceDictionary R= E.coli_K12_MG1655.fa O= E.coli_K12_MG1655.dict
得到变异文件
gatk HaplotypeCaller -R E.coli_K12_MG1655.fa -I SRR1770413.dup.bam -O SRR1770413.dup.vcf
查看变异位点
samtools tview SRR1770413.dup.bam E.coli_K12_MG1655.fa
筛选vcf，得到变异位点
vcftools --vcf SRR1770413.dup.vcf --max-missing 0.5 --minDP 3 --minQ 30 --recode --recode-INFO-all --out filt
对vcf文件进行划窗处理
vcftools --vcf filt.recode.vcf --window-pi 100 --window-pi-step 5
sed 's/\t/,/g' out.windowed.pi> a.csv

其中使用到了以下软件：
SRA Tookit
bgzip：用于解压缩gzip格式的压缩文件。
BWA ：用于建立参考基因组的索引以及进行序列比对。
samtools：一套用于处理SAM格式文件的工具，这里用于建立FAI索引和质量检查。
picard：Picard是一个用于处理BAM文件的工具集。
gatk ：用于基因组数据分析的工具包。
vcftools：用于处理VCF文件的工具。
建议在环境中安装conda以更加简便地管理软件。