WGS (whole genome sequencing) 全基因组测序分析，以Saccharomyces cerevisiae SY14菌株的illumin 测序数据为例。
1. rawdata数据处理
下载数据
SY14: SRR6825081； BY4742 (control): SRR6825082
以SRR6825081为例，这里软件都调用的绝对路径。
#创建脚本download_data.sh
vim download_data.sh
#!/bin/bash
sratoolkit="/sibcb/program/install/sratoolkit.2.10.0-centos_linux64/bin/"
$sratoolkit/prefetch-orig.2.10.0 `$sratoolkit/srapath-orig.2.10.0 SRR6825081`
#提交脚本
qsub -N download -cwd download_data.sh
sra格式转fastq格式
sra格式是二进制文件，Sequence Read Archive (SRA) 数据库，储存来自不同高通量测序平台的原始序列数据和比对信息。
#创建脚本sra2fastq.sh
vim sra2fastq.sh
#!/bin/bash
/sibcb/program/install/sratoolkit.2.10.0-centos_linux64/bin/fastq-dump-orig.2.10.0 -O ~/WGS/SY14/illumina/fastq --split-3 --gzip ~/WGS/SY14/illumina/SRA/SRR6825081/SRR6825081.sra
#提交脚本sra2fastq.sh
qsub -N fastq -cwd sra2fastq.sh
双端数据生成两个文件：SRR6825081_1.fastq.gz，SRR6825081_2.fastq.gz
利用fastqc对测序数据进行质量评估（Quality Control）
#创建脚本fastqc.sh
vim fastqc.sh
#!/bin/bash
/sibcb/program/install/FastQC-11.8/fastqc ~/WGS/SY14/illumina/fastq/SRR6825081_1.fastq.gz ~/WGS/SY14/illumina/fastq/SRR6825081_2.fastq.gz -o fastq
#提交脚本fastqc.sh
qsub -N fastqc -cwd fastqc.sh
trim_galore 过滤adapter, cut低质量的碱基
#创建脚本 trim_galore.sh
vim trim_galore.sh
#!bin/bash
sibcb/program/install/TrimGalore-0.5.0/trim_galore --output_dir ~/WGS/SY14/illumina/clean_data -q 25 --phred33 --length 36 --paired ~/WGS/SY14/illumina/fastq/SRR6825081_1.fastq.gz ~/WGS/SY14/illumina/fastq/SRR6825081_2.fastq.gz
#提交脚本 trim_galore.sh
qsub -N trim_galore -cwd trim_galore.sh -q g5.q
2. 序列比对、排序
ssh，随机分配一个服务器节点后，运行如下命令行，logout退出，返回到登录节点。
bwa index 构建索引
/sibcb/program/install/bwa-0.7.15/bwa index ~/rnaseq/1_rnaseq/yeast_genome/GCF_000146045.2_R64_genomic.fna
bwa mem比对
/sibcb/program/install/bwa-0.7.15/bwa mem -M -t 10 -R '@RG\tID:SRR6825081\tSM:SRR6825081\tLB:NGS\tPL:Illumina' GCF_000146045.2_R64_genomic.fna SRR6825081_1_val_1.fq.gz SRR6825081_2_val_2.fq.gz > SY14_illumina.bam 2 > SY14.mem.log
-R 设置reads标头，ID：通道名或样本名,通过这个信息分组，必须唯一；SM：样本名；LB：文库名；PL：测序平台信息
">" 标准输出到SY14_illumina.bam；"2 >" 标准误输出到SY14.mem.log
samtools排序
samtools sort SY14_illumina.bam -o SY14_illumina.bam.sort
#排序
samtools index SY14_illumina.bam.sort
#构建索引
3. 利用GATK（Genome Analysis Toolkit）call SNP、INDEL
GATK是用java语言写的，调用时注意java版本。
gatk MarkDuplicates 标记重复序列
/sibcb/program/install/gatk-4/gatk MarkDuplicates -I SY14_illumina.bam.sort -O SY14_bam.sort.markdup.bam -M SY14_bam.sort.markdup_metrics.txt
#标记重复序列
samtools index SY14_bam.sort.markdup.bam
#生成.bai索引文件
gatk CreateSequenceDictonary 生成字典文件.dict
/sibcb/program/install/gatk-4/gatk CreateSequenceDictionary -R GCF_000146045.2_R64_genomic.fna -O GCF_000146045.2_R64_genomic.dict
#生成字典文件.dict
samtools faidx GCF_000146045.2_R64_genomic.fna
#为fna文件生成fai文件
gatk HaplotypeCaller 输出vcf
/sibcb/program/install/gatk-4/gatk HaplotypeCaller -R GCF_000146045.2_R64_genomic.fna -I SY14_bam.sort.markdup.bam -O SY14.vcf
生成的vcf表格不同列，详细讲解可参考：
Input & output files - SnpEff & SnpSift Documentation (pcingola.github.io)
4. snpEff 注释
snpEff的注释，在vcf文件的第八列后面加入 ANN 详细信息。
snpEff 运行可参考：
Running SnpEff - SnpEff & SnpSift Documentation (pcingola.github.io)
cd snpEff
mkdir data
cd data #在data文件夹下分别创建genomes文件夹、Sc文件夹
mkdir genomes
mkdir Sc
cd genomes
cp ~/WGS/SY14/illumina/clean_data/GCF_000146045.2_R64_genomic.fna Sc.fa
#重命名
cd Sc
cp ~/WGS/SY14/illumina/clean_data/GCF_000146045.2_R64_genomic.gff genes.gff
cd ~/WGS/SY14/illumina/clean_data/anotation/snpEff
echo "Sc.genome:Sc">>snpEff.config
#配置snpEff.conf文件
java -jar snpEff.jar build -c snpEff.config -gff3 -v Sc
#手动配置数据库
java -jar: Java环境下运行程序
-c：snpEff.config配置文件路径
-gff3：设置输入基因组注释信息是gff3格式
-v：设置在程序运行过程中输出的日志信息
Sc参数：即设置输入的基因组版本信息。
cd ~/WGS/SY14/illumina/clean_data/anotation/snpEff
java -jar snpEff.jar -v Sc ~/WGS/SY14/illumina/clean_data/SY14/SY14.vcf > SY14.annotation.vcf
Sc表示调用 data/Sc/中的注释信息；
最终生成三个文件，SY14.annotation.vcf（查看结果）；snpEff_genes.txt；snpEff_summary.html（类似fastqc网页版的总结）
网友评论