使用gmap 进行基因组结构预测

2023-07-13  本文已影响0人  惊鸿影

可以直接用conda安装

conda install -c bioconda gmap

如下步骤假设你有一个基因组序列和一个参考的同种或近缘CDS序列,分别命名为"genome.fa"和"cds.fa" (其中genome是作为reference来构建索引的,和cds一定要区分开,前面看到的好多教程都没说清哪个作为reference的,弄混了好久)

第一步:构建GMAP/GSNAP索引数据库

GMAP对FASTA文件中每个记录下的序列的长度有一定限制, 每一条不能超过4G, 能应付的了大部分物种了。

构建索引分为两种情况考虑,第一种是一个fasta文件包含所有的序列

gmap_build -D ./genome -d genome_reference genome.fa

# -D 指定建立索引的文件夹的位置,可以根据你的基因组名字新建一个目录
# -d 指定建立索引的文件夹的名字  在这里我是新建了一个genome的目录再在genome的目录下新建了一个genome_reference目录作为存放索引

第二种则是每个染色体的序列都单独存放在一个文件夹里,比如说你下载人类参考基因组序列解压后发现有N多个fasta文件, 然后你就想用其中几条染色体构建索引

gmap_build  -D ./genome -d genome_reference Chr1.fa 
Chr2.fa Chr3.fa ...

注: 这里的-d表示数据K库的名字,默认把索引存放在gmap安装路径下的share里,可以用-D更改.此外还有一个参数-k用于设置K-mer的长度, 默认是15, 理论上只有大于4GB基因组才会有两条一摸一样的15bp序列(当然是完全随机情况下)。

gmap -D ./genome_gmap -d genome_reference   -t 8 -f  gff3_gene cds.fa > genome.gff3
#千万不要用nohup运行

之后获得一个gff3文件,就可以使用gffread提取CDS区域了

根据GFF或者GTF提取CDS,蛋白质和外显子序列

gffread genome.gff3 -g genome.fa -x cds.fa -y pep.fa -w cdna.fa

只提取翻译后蛋白序列

gffread genome.gff3 -g genome.fa -y tr_pep.fa

根据reference提取CDS序列

gffread genome.gff3 -g genome.fa -x cds.fa

只提取外显子序列

gffread genome.gff3 -g genome.fa -w exons.fa
上一篇下一篇

猜你喜欢

热点阅读