基因组数据准备
在进行数据分析之前，需要准备好参考基因组数据，基因组数据可以从公共数据库下载，也可以根据基因组文献提供的地址到指定网站进行下载。一套完整的数据至少包括如下内容：
- 基因组序列文件，fasta 格式
- 基因结构注释文件，gtf/gff3 格式
- 蛋白质序列文件，fasta 格式
- cds 序列文件，fasta 格式
BSA 分析使用的参考基因组需要组装到染色体水平，下面以拟南芥基因组为例进行演示。
拟南芥基因组下载自 Ensembl Plants 数据库。
数据库网址：http://plants.ensembl.org/index.html。
# 数据下载
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-56/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.56.gff3.gz
# 解压缩
gzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.56.gff3.gz
数据统计
# 统计基因组整体信息
seqkit stats Arabidopsis_thaliana.TAIR10.dna.toplevel.fa > genome.allstat
# 统计每条序列长度
seqkit fx2tab -g -l -n -i Arabidopsis_thaliana.TAIR10.dna.toplevel.fa > genome.seqstat
# 统计编码基因信息
awk -F "\t" '$3=="gene" && $NF ~ \
/protein_coding/{ num+=1; total += $5-$4+1;}END{print "number: \
"num"\ntotal_len: "total "\navg_len: "total/num}' Arabidopsis_thaliana.TAIR10.56.gff3 > gene.stat
# 统计编码转录本信息
awk -F "\t" '$3=="mRNA" && \
$NF ~ /protein_coding/{ num+=1; total += $5-$4+1;}END{print "number: \
"num"\ntotal_len: "total "\navg_len: "total/num}' \
Arabidopsis_thaliana.TAIR10.56.gff3 > transcript.stat
结果如下
$ cat genome.allstat
file format type num_seqs sum_len min_len avg_len max_len
genome.fa FASTA DNA 7 119,667,750 154,478 17,095,392.9 30,427,671
$ cat genome.seqstat
1 30427671 35.68
2 19698289 35.86
3 23459830 36.32
4 18585056 36.20
5 26975502 35.93
Mt 366924 44.77
Pt 154478 36.29
$ cat gene.stat
number: 27628
total_len: 65554418
avg_len: 2372.75
$ cat transcript.stat
number: 48321
total_len: 128224748
avg_len: 2653.6
功能注释
我们使用在线版 eggnog-mapper 进行 eggnog 数据库注释，在线版 eggnog-mapper 单次输入序列条数不能超 10 万条，如果序列条数超过 10 万条，需要进行切分。蛋白数据准备好以后，登录 http://eggnog-mapper.embl.de/ 上传蛋白序列, 在线进行蛋白注释，完成后下载注释结果 *emapper.annotations。
网友评论