一、下载数据

1. 下载参考基因组

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/350/185/GCA_000350185.1_ASM35018v1/GCA_000350185.1_ASM35018v1_genomic.fna.gz

gunzip GCA_000350185.1_ASM35018v1_genomic.fna.gz

2. 下载SAR数据

Wget https://sra-pub-run-odp.s3.amazonaws.com/sra/DRR398234/DRR398234 #下载SAR数据集

因为下载下来的文件并非输入文件格式fasta格式，且此sar数据为双端测序，若不拆分开，则在后续的比对中易出现问题，所以还要对下载下来的sar数据进一步处理。

mv DRR398234 DRR398234.sra

fastq-dump

fastq-dump --split-3 DRR398234.sra

结果生成DRR398234_1.fastq DRR398234_2.fastq

二、使用BWA进行数据比对

bwa index GCA_000350185.1_ASM35018v1_genomic.fna -p 398200_index #为参考基因组建立索引

bwa mem 398200_index DRR398234_1.fastq DRR398234_2.fastq>test_bwa_398200.sam  #进行bwa比对

less test_bwa_398200.sam #查看

图片.png

使用BWA进行比对后可以得到一个sam文件，是以TAB为分割符的文本格式。

sam文件解读：
@SQ，参考序列说明；
@PG，使用的程序说明；
“DRR398234.1”：read序列的名字；
“83”=64+16+1+2：概括出一个合适的标记，(1：序列是一对序列中的一个，2：比对结果是一个pair-end比对的末端，16：在这个比对上的位点，序列与参考序列反向互补，64：序列是mate1)；
“AP012306.1”：参考序列的名字；
“1828021”：在参考序列上的位置；
“60”：MAPQ值，越高则越独特，此处数值较高，则独特性较高；
“296M”：代表比对结果的CIGAR字符串，这段字符的意思是296个匹配；
“=”：成对的另一个片段比对上的参考序列的编号，此处为=，则为同一个片段，若没有另外的片段，则是；
“1828021”：成对的另一个片段在参考序列上的位置；如果不可用，此处为0；
“-296”：估计出的片段的长度，此处为负值，则mate序列位于本序列的上游，不可用时为0；
“TCAAAT...”：为SEQ，read的序列；
“6544)<4<::<92D>=....”：QUAL，为ASCll码格式的序列质量；
后续参数暂且不做详细解释。
此为第一条比对结果，后续结果可以参照此条进行解读。

三、使用samtools转换文件格式，查看比对结果

samtools faidx GCA_000350185.1_ASM35018v1_genomic.fna #为参考基因组建立索引，生成了prefix.fai文件

less GCA_000350185.1_ASM35018v1_genomic.fna.fai #查看

图片.png

第一列为染色体编号，第二列为染色体长度，第三列为所占的字符数，第四列为每一行的碱基个数，第五列为每一行的占位符。

samtools view -bhS -t GCA_000350185.1_ASM35018v1_genomic.fna.fai -o 398200_bwa.bam test_bwa_398200.sam #sam文件转为bam文件

less 398200_bwa.bam #查看

图片.png

sam文件和bam文件都可以存储序列比对信息，但是bam是以二进制存储，可以节约空间，且计算机可读。我对将sam文件转为bam文件没有什么要求，所以我直接沿用了老师教学时使用的参数，正因为bam文件为二进制文件，所以看上去像是一堆乱码，其实不是，但是没关系，我们不需要去解读它，也不会影响后续操作。

samtools sort 398200_bwa.bam -o 398200_bwa.bam.sorted #为bam文件排序，sort只能为bam文件排序，而不能为sam；不同版本samtools sort命令的-o参数不同

samtools index 398200_bwa.bam.sorted #为排序的bam文件建立索引. *.bai文件

samtools tview 398200_bwa.bam.sorted GCA_000350185.1_ASM35018v1_genomic.fna #查看比对结果

此处由于我对使用samtools不熟练，所以跳过此步，不进行查看了，直接进行下一步。

samtools depth 398200_bwa.bam.sorted>398200_depth.txt #测试参考基因组每个位点或一段区域的测序深度

less 398200_depth.txt

图片.png

从左到右依次为参考基因组序列名字、在基因组上的位置、深度值

samtools flagstat 398200_bwa.bam.sorted #统计比对结果

图片.png

解读：总共有2021505个reads，总体上reads的匹配率为99.92%（匹配率很高），有2018572个reads属于paird reads,reads1中有1009286个reads，reads2中有1009286个reads,有2008940个reads正确配对，一对reads竣比对上的reads数为2016078，只有单条reads比对上的reads数为850。

samtools mpileup -f GCA_000350185.1_ASM35018v1_genomic.fna 398200_bwa.bam.sorted >bcf_2.txt

less bcf_2.txt

从此文件我们可以得知参考序列名、位置、参考碱基、比对上的reads数、比对情况以及比对上碱基的质量。

图片.png

samtools mpileup -gf GCA_000350185.1_ASM35018v1_genomic.fna 398200_bwa.bam.sorted>398200_bwa.bcf #生成bcf文件

less 398200_bwa.bcf

bcf文件也为二进制文件

图片.png

四、使用bcftolls检测变异位点

bcftools call -vm 398200_bwa.bcf -o 398200_bwa.variant.bcf

bcftools view -v snps,indels 398200_bwa.variant.bcf > 398200_bwa.snps.vcf

less 398200_bwa.snps.vcf

图片.png

此时我们可以得到参考序列该位置的allele，variant的信息等，

bcftools filter -o 398200_bwa.snps.filtered.vcf -i 'QUAL>20 && DP>5' 398200_bwa.snps.vcf #变异位点过滤

less 398200_bwa.snps.filtered.vcf

对于变异位点的过滤，我依旧沿用了老师的命令及参数。

图片.png

五、使用Annovar注释变异位点

convert2annovar.pl -format vcf4 398200_bwa.snps.filtered.vcf > 398200_bwa.snps.filtered.avinput #生成annovar输入文件

less 398200_bwa.snps.filtered.avinput

图片.png

从左到右依次为参考基因序列名字，起始位置，末端位置，参考核苷酸，观察到的核苷酸，其他的我们暂且不需要注意。

/disk1/shares/gff3ToGenePred

/disk1/shares/gff3ToGenePred -useName GCA_000350185.1_ASM35018v1_genomic.gff 3982-genome_refGene.txt

less 3982-genome_refGene.txt

图片.png

cut-f 12 3982-genome_refGene.txt > column_3982.txt

cut-f 2-15 3982-genome_refGene.txt > column_3982else.txt

paste column_3982.txt column_3982else.txt >3982-genome-new_refGene.txt

retrieve_seq_from_fasta.pl -format refGene -seqfile GCA_000350185.1_ASM35018v1_genomic.fna-outfile 3982-genome-new_refGeneMrna.fa 3982-genome-new_refGene.txt

less 3982-genome-new_refGeneMrna.fa

图片.png

cp 3982-genome-new_refGene* ~/Biosofts/annovar/humandb/

ll ~/Biosofts/annovar/humandb/  #自定义注释数据库

由于我是自己在网上寻找的实验数据，所以最好是自己定义一个注释数据库，主要需要用到的就是参考基因组的gff文件，gff文件主要就是用来注释基因组的，代码参考老师在教学时使用的。

annotate_variation.pl --geneanno --dbtype refGene --buildver 3982-genome-new 398200_bwa.snps.filtered.avinput ~/Biosofts/annovar/humandb/ #注释变异位点

less 398200_bwa.snps.filtered.avinput.exonic_variant_function

图片.png

此为avinput.exonic_variant_function文件，第一列为该变异在输入文件的行号，第二列为对编码基因的影响，第三列为被影响的基因或转录本，其余列则与输入文件相同。

less 398200_bwa.snps.filtered.avinput.variant_function

图片.png

此为avinput.variant_function文件，其中第一列为变异分类，第二列为基因名或Symbol，后面列则与avinput输入文件相同。