可以直接用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
网友评论