基因组数据的重测序分析
我通过查资料获得已知达松维尔拟诺卡氏菌亚种（cardiopsis dassonvillei subsp. dassonvillei）的基因组原始测序序列，我想知道这个亚种与达松维尔拟诺卡氏菌（Nocardiopsis dassonvillei ）的基因组相比有哪些不同。
1.需要的软件
软件名：Aspera 版本号：3.6.2.117442
软件名：sratoolkit 版本号：2.9.2
软件名：FastQC 版本号：0.11.7
软件名：Trimmomatic版本号：0.38
软件名：bwa 版本号：0.7.17-r1188
软件名：samtools 版本号：1.7
软件名：Annovar 版本：$Date: 2017-07-17 01:16:48 -0400 (Mon, 17 Jul 2017) $
2.数据下载
达松维尔拟诺卡氏菌亚种（cardiopsis dassonvillei subsp. dassonvillei）的基因组：
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR022/SRR022534/SRR022534.sra
~/.aspera/connect/bin/ascp -T -i /home/lizeguo/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR022/SRR022534/SRR022534.sra ./Seqs/
达松维尔拟诺卡氏菌（Nocardiopsis dassonvillei ）参考基因组序列：
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz
~/.aspera/connect/bin/ascp -T -i /home/lizeguo/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m anonftp@ftp-private.ncbi.nlm.nih.gov:/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz ~/Seqs/
基因组gff文件
cd ~/Seqs
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.gff.gz
3.主要分析步骤和结果
（一）序列与参考基因组的比对
1.数据文件的格式转换
在NCBI的SRA数据库检索序列号SRR022534查看测序数据类型为454的单末端测序
2018-12-20 10-57-15屏幕截图.png
fastq-dump SRR022534.sra
2018-12-20 11-04-11屏幕截图.png
2.质量评估
mkdir fastqc_result
fastqc SRR022534.fastq
mv *.zip *.html fastqc_result/
2018-12-20 11-16-21屏幕截图.png
3.测序数据的数据过滤
mkdir trim_out
java -jar ~/BioSofts/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 SRR022534.fastq ./trim_out/SRR011534_out.fastq.gz ILLUMINACLIP:/home/lizeguo/BioSofts/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75
2018-12-20 11-29-33屏幕截图.png
4.建立参考基因组索引
gunzip GCA_001877055.1_ASM187705v1_genomic.fna.gz
bwa index GCA_001877055.1_ASM187705v1_genomic.fna
2018-12-20 11-37-59屏幕截图.png
5.测序数据比对到参考基因组得到sam文件
bwa mem GCA_001877055.1_ASM187705v1_genomic.fna trim_out/SRR011534_out.fastq.gz >bwa_mem_SRR011534.sam
6.sam文件转换为bam文件
samtools faidx GCA_001877055.1_ASM187705v1_genomic.fna
samtools view -bhS -t GCA_001877055.1_ASM187705v1_genomic.fna.fai -o bwa_mem_SRR011534.bam bwa_mem_SRR011534.sam
7.为bam文件排序
samtools sort bwa_mem_SRR011534.bam -o bwa_mem_SRR011534.sorted.bam
8.为bam文件建立索引
samtools index bwa_mem_SRR011534.sorted.bam
9.显示基因组比对情况
samtools tview bwa_mem_SRR011534.sorted.bam GCA_001877055.1_ASM187705v1_genomic.fna
reads比对情况
10.测试参考基因组每个位点或一段区域的测序深度
samtools depth bwa_mem_SRR011534.sorted.bam >>depth.txt
less depth.txt
每个位点或区域的测序深度
11.统计比对结果
samtools flagstat bwa_mem_SRR011534.sorted.bam
比对结果
结果显示有458257个碱基，有441750个碱基匹配上了，占比96.40%.
（二）变异位点的检测
1.去除PCR重复
samtools rmdup bwa_mem_SRR011534.sorted.bam bwa_mem_SRR011534_nopcr.bam
2.生成bcf文件
samtools mpileup -gf GCA_001877055.1_ASM187705v1_genomic.fna bwa_mem_SRR011534_nopcr.bam >bwa_mem_SRR011534.bcf
3.基因变异检测，得到bwa_mem_SRR011534.snps.bcf文件
bcftools call -vm bwa_mem_SRR011534.bcf -o bwa_mem_SRR011534.variants.bcf
bcftools view -v snps,indels bwa_mem_SRR011534.variants.bcf > bwa_mem_SRR011534.snps.vcf
less bwa_mem_SRR011534.snps.vcf
2018-12-31 13-47-46屏幕截图.png
4.变异位点的过滤
bcftools filter -o bwa_mem_SRR011534.snps.filtered.vcf -i 'QUAL>20 &&DP>5' bwa_mem_SRR011534.snps.vcf
（三）变异基因注释
1.生成annovar输入文件
convert2annovar.pl -format vcf4 bwa_mem_SRR011534.snps.vcf > bwa_mem_SRR011534.snps.avinput
2.自定义注释数据库
gunzip GCA_001877055.1_ASM187705v1_genomic.gff.gz
2.1gff3文件转为GenePred文件
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/gff3ToGenePred
chmod 777 gff3ToGenePred
./gff3ToGenePred -useName GCA_001877055.1_ASM187705v1_genomic.gff 7055-genome_refGene.txt
2.2修改GenePred文件
cut -f 12 7055-genome_refGene.txt >column1.txt
cut -f 2-15 7055-genome_refGene.txt >column_else.txt
paste column1.txt column_else.txt >7055-genome_new_refGene.txt
2.3为建立每个基因与编码序列对应文件
retrieve_seq_from_fasta.pl -format refGene -seqfile GCA_001877055.1_ASM187705v1_genomic.fna -outfile 7055-genome_new_refGeneMrna.fa 7055-genome_new_refGene.txt
2.4拷贝数据库文件到annovar安装目录humandb文件夹
cp 7055-genome_new_refGene* ~/BioSofts/annovar/humandb/
3.注释变异基因位点，生成avinput.variant_function和avinput.exonic_variant_function后缀的两个结果文件
annotate_variation.pl --geneanno --dbtype refGene --buildver 7055-genome_new bwa_mem_SRR011534.snps.avinput ~/BioSofts/annovar/humandb/
2018-12-31 15-13-41屏幕截图.png
4.查看结果
less bwa_mem_SRR011534.snps.avinput.variant_function
2018-12-31 15-15-49屏幕截图.png
less bwa_mem_SRR011534.snps.avinput.exonic_variant_function
2018-12-31 15-15-49屏幕截图.png
