1、软件与环境

（1）环境
    #操作系统
        ubuntu、deepin
    #java
        #试了一下，bcftools和fastq在java1.8上运行不了
        #所以改成默认的java版本
        #sudo apt-get install openjdk-8-jdk
        sudo apt-get install default-jre
        #如果安装了多个java版本，以下命令切换成java-10
        update-alternatives --config java

（2）软件
    sratoolkit
        fastq-dump #解压文件
        prefetch #下载ncbi sra文件
    fastqc
    bwa 
    samtools
    trimmomatic
    bcftools
    #安装bwa samtools fastaqc bcftools
    sudo apt install bwa samtools fastaqc bcftools
    #stratoolkit和trimmoatic引用命令所在位置来执行
(3)环境变量配置

2、数据

SRR1770413_1.sra                      #测序文件
GCF_000005845.2_ASM584v2_genomic.fna  #参考序列

3、步骤

1、#解压sra文件：fastq-dump
soft/sratoolkit.2.8.2-ubuntu64/bin/fastq-dump  --split-files SRR1770413_1.sra

2、#fastqc质量控制：fastqc
mkdir fastqc_result
fastqc SRR1770413_1*.fastq -o fastqc_result

#除去接头和低质量序列：Trimmomatic
  #压缩成gz文件
  gzip SRR1770413_1_1.fastq
  gzip SRR1770413_1_2.fastq
  #去除接头序列
  java -jar soft/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 -trimlog logfile SRR1770413_1_1.fastq.gz SRR1770413_1_2.fastq.gz -baseout filtered.fq.gz ILLUMINACLIP:soft/Trimmomatic-0.38/adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50

3、#一致性序列的获得
  #用bwa将序列比对回基因组，因为是双端数据，比对完后还要合并一下
  #解压文件
  gzip -d GCF_000005845.2_ASM584v2_genomic.fna.gz
  gunzip filtered_1P.fq.gz filtered_2P.fq.gz 
  #建立索引
  bwa index GCF_000005845.2_ASM584v2_genomic.fna
  #采用8线程比对回基因组
  bwa aln -t 8 GCF_000005845.2_ASM584v2_genomic.fna filtered_1P.fq > read1.sai
  bwa aln -t 8 GCF_000005845.2_ASM584v2_genomic.fna filtered_2P.fq > read2.sai
  #合并比对结果为sam文件
   bwa sampe GCF_000005845.2_ASM584v2_genomic.fna read1.sai read2.sai filtered_1P.fq filtered_2P.fq > reads.sam
  #将sam格式转化bam格式
  samtools view -bS reads.sam > reads.bam
  #排序 
  samtools sort reads.bam > reads.sort.bam

  #将bam转为bcf格式
  bcftools mpileup -f GCF_000005845.2_ASM584v2_genomic.fna reads.sort.bam > reads.bcf
  #将bcf转为vcf格式
  bcftools call -c reads.bcf > reads.vcf
  #将vcf压缩成gz格式，有多种方式
  bcftools view reads.seq.vcf -Oz -o reads.seq.vcf.gz
  #建立索引
  bcftools index reads.vcf.gz
  #得到一致性序列
  bcftools consensus -f GCF_000005845.2_ASM584v2_genomic.fna reads.vcf.gz > reads.fa
4、snp、indel分析
  #使用的是samtools的mpileup
  #bcftools mpileup -t DP,SP -ugf GCF_000005845.2_ASM584v2_genomic.fna reads.sort.bam > reads_snp.bcf
  samtools mpileup -t DP,SP -ugf GCF_000005845.2_ASM584v2_genomic.fna reads.sort.bam > reads_snp.bcf
  bcftools call -vm reads.bcf > reads.variant.vcf
  bcftools call -V indels -vm reads.bcf > reaads.snps.vcf
  bcftools call -V snps -vm reads.bcf > reaads.indels.vcf