基因组测序分析linux,基因组数据的重测序分析

我通过查资料获得已知达松维尔拟诺卡氏菌亚种(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.数据下载

~/.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/

~/.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