此次实验操作由于参考基因组下载出问题，实验中止，正确操作请参考下一篇。

一、下载数据

1、下载SAR数据

进入NCBI
whole genome sequencing - BioProject - NCBI (nih.gov)

图片.png

选择[Escherichia coli]

点击SRA Experiments后的8

随便挑选一条

点击第一个*All runs*

图片.png

1.

图片.png

2.点击下面表格中的SRR号

234

图片.png

图片.png

切换到‘Data access’界面，就找到数据链接了

图片.png

复制链接

用wget下载

wget https://sra-pub-run-odp.s3.amazonaws.com/sra/DRR398233/DRR398233

下载

图片.png

2.下载参考基因

图片.png

可以在STRAIN找到它的菌株种型，然后复制名称，进入NCBI中搜索

图片.png

如图，第一个结果为最保险的参考基因，点击进入

图片.png

图片.png

图片.png

图片.png

为参考基因组建立索引

bwa index sequence.fasta -p 398200_index

图片.png

二、bwa比对

bwa mem 398200_index  DRR398233 DRR398234 >test_bwa_398200.sam

图片.png

less test_bwa_398200.sam

图片.png

 mv DRR398234 DRR398234.sra

图片.png

fastq-dump
fastq-dump DRR398234.sra

图片.png

将文件格式从sar转为fastq

图片.png

ll #查看文件夹内容

图片.png

bwa index sequence.fasta -p 398200_index #为参考基因组建立索引
bwa mem 398200_index DRR398234.fastq >test_bwa_398200.sam #bwa比对

图片.png

less test_bwa_398200.sam #查看

图片.png

fastq-dump --split-3 DRR398234.sra

图片.png

结果生成DRR398234_1.fastq
DRR398234_2.fastq

再次比对

bwa mem 398200_index DRR398234_1.fastq DRR398234_2.fastq>test_bwa_398200.sam

图片.png

less test_bwa_398200.sam

图片.png

三、 samtools

SAM格式转为BAM格式

samtools faidx sequence.fasta
  #为参考基因组建立索引，生成了prefix.fai文件

less sequence.fasta.fai

按q退出

samtools view -bhS -t sequence.fasta.fai -o 398200_bwa.bam test_bwa_398200.sam #(最后这个sam文件用大家自己目录的sam文件名)   ##sam文件转为bam文件

图片.png

less 398200_bwa.bam

按y

按q退出

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

图片.png

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

图片.png

统计比对结果

samtools depth 398200_bwa.bam.sorted>398200_depth.txt

图片.png

less 398200_depth.txt

按q退出

samtools flagstat 398200_bwa.bam.sorted

图片.png

查看比对结果

samtools tview

samtools tview 398200_bwa.bam.sorted sequence.fasta

图片.png

g, 输入：CP000100.1:1000

此处由于我不知道输入应为什么信息，所以无法呈现......所以跳过此步

BAM2BCF

samtools mpileup -f sequence.fasta 398200_bwa.bam.sorted >bcf_2.txt

图片.png

less bcf_2.txt

图片.png

samtools mpileup -gf sequence.fasta 398200_bwa.bam.sorted>398200_bwa.bcf

图片.png

less 398200_bwa.bcf

按y

图片.png

bcftools

bcftolls检测变异位点

cd ~/bwa_test

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

图片.png

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

图片.png

less 398200_bwa.snps.vcf

按q退出

变异位点过滤

bcftools filter -o 398200_bwa.snps.filtered.vcf -i 'QUAL>20 && DP>5' 398200_bwa.snps.vcf

bcftools filter [OPTIONS] FILE 该命令用于固定阈值过滤

-i, --include EXPRESSION 接表达式命令，当时表达式返回为真值时，输出该位点。表达式命令和使用见表达式部分

[最全bcftools使用说明--只看本文就够了 - 简书 (jianshu.com)]https://www.jianshu.com/p/99895c7338b2

图片.png

less 398200_bwa.snps.filtered.vcf

图片.png

五、使用Annovar注释变异位点

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

图片.png

less 398200_bwa.snps.filtered.avinput

图片.png

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/008/248/145/GCA_008248145.1_ASM824814v1/GCA_008248145.1_ASM824814v1_genomic.gbff.gz  #下载gbff文件

图片.png

gunzip GCA_008248145.1_ASM824814v1_genomic.gbff.gz #解压

图片.png

perl bp_genbank2gff3.pl GCA_008248145.1_ASM824814v1_genomic.gbff  #将gbff格式转为gff格式

失败，没有bp_genbank2gff3.pl文件，在网上搜索，下载了一个
https://blog.csdn.net/weixin_42677653/article/details/106977642

less GCA_008248145.1_ASM824814v1_genomic.gbff.gff

图片.png

再次尝试

perl bp_genbank2gff3.pl GCA_008248145.1_ASM824814v1_genomic.gbff

图片.png

图片.png

失败，再次重新找

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

图片.png

图片.png

gunzip GCA_000350185.1_ASM35018v1_genomic.gff.gz  #解压

图片.png

/disk1/shares/gff3ToGenePred

图片.png

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

图片.png

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

图片.png

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

图片.png

retrieve_seq_from_fasta.pl -format refGene -seqfile sequence.fasta -outfile 3982-genome-new_refGeneMrna.fa 3982-genome-new_refGene.txt

图片.png

图片.png

less 3982-genome-new_refGeneMrna.fa

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

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

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

less 398200_bwa.snps.filtered.avinput.variant_function