根据已知的参考基因组和GTF创建适合自己分析的基因组和GTF
###############################################################
######################## Genome build #########################
###############################################################
###删除染色体号后面的内容(包括空格)
###写了一个perl程序,可以切除匹配项之后的信息,得到匹配上的项。
###Perl 名 :cut_space.pl,以下就是程序内的语句:(不同文件需要修改染色体号)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | #!/usr/bin/perl open(FR,"10.fa")||die; open(FW,">10_cut.fa")||die; while(<FR>) { chomp; if ($_=~/(>\d+)/) { print FW "$1\n"; } else {print FW "$_\n";} } close FR; close FW; |
###有些时候染色体号前需要加上chr,可以用perl程序:
###程序名:fa_add_chr.pl:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | #!/usr/bin/perl open(FR,"10.fa")||die; open(FW,">10_add_chr.fa")||die; while(<FR>) { chomp; if($_ =~/^>/) { $_ =~ s/>/>chr/g; print FW "$_\n" } else {print FW "$_\n";} } close FR; close FW; |
###############################################################
######################## Ref GTF build ########################
###############################################################
###两种方法:一种是直接通过perl实现:
###一、perl程序:两步;
###1、删除没用的染色体:gtf_del_useless_chr.pl
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | #!/usr/bin/perl open(FW,"Macaca_mulatta.MMUL_1.73.gtf")||die; open(FM,">DEL_Macaca_mulatta.MMUL_1.73_del_un.gtf")||die; while(<FW>) { chomp; if($_=~/^1099*/) { print FM "" } else { print FM "$_\n" } } close FW; close FM; |
###2、在染色体号前加chr:gtf_add_chr.pl
1 2 3 4 5 6 7 8 9 10 11 12 13 | #!/usr/bin/perl #first chr numble must add chr. open(FR,"DEL_Macaca_mulatta.MMUL_1.73_del_un.gtf")||die; open(FW,">Add_chr_Macaca_mulatta.MMUL_1.73_del_un.gtf")||die; while(<FR>) { chomp; if($_=~/^\d/ or $_=~/^X/ or $_=~/^MT/) {print FW "chr$_\n"} else {print FW "\n"} } close FR; close FW; |
###二、Linux shell实现:
###1、导出每个染色体的GTF:
1 2 3 4 5 6 7 8 9 | for i in `seq 20` do grep '^'${i}''$'\t' Macaca_mulatta.MMUL_1.73.gtf >${i}.gtf done grep '^X'$'\t' Macaca_mulatta.MMUL_1.73.gtf >X.gtf grep '^MT'$'\t' Macaca_mulatta.MMUL_1.73.gtf >MT.gtf cat *.gtf >Macaca_mulatta_1.73_final.gtf #合并每个GTF; |
###2、在染色体号前加chr:gtf_add_chr.pl,同perl的方法相同。(略)
相关推荐:
- NCBI在线BLAST使用方法与结果详解 2938
- 神经网络术语:Epoch、Batch Size和迭代 527
- Consed的安装与使用教程 465
- 陈连福的NGS生物信息学培训教材V2.1 277
- WGCNA分析使用教程 272
-
***3zero_jason 2014-09-11 10:33#2
我的参考基因组名称为,indica_dna_ref,我首先使用bowtie2-build建立index。我的命令是 tophat -p 4 -G indica_ref.gff -o DGDZ_OUT indica_ref diguduizhao_1.fq diguduizhao_2.fq.tophat结果的log中的g2f.err的说明是:GFF Error:duplicate/invalid ‘transcrip’ feature ID=transcript:BGIOSGA002656-TA. 另外,我已经把基因组和GFF文件格式整理成第一列为Chr*的格式了。[img]http://[/img]
-
***4baker 2015-08-20 22:10#3
你是怎么解决的啊,我也遇到这个问题
最新创建圈子
-
原料药研发及国内外注册申报
2019-01-25 10:41圈主:caolianhui 帖子:33 -
制药工程交流
2019-01-25 10:40圈主:polysciences 帖子:30 -
健康管理
2019-01-25 10:40圈主:neuromics 帖子:20 -
发酵技术
2019-01-25 10:39圈主:fitzgerald 帖子:17 -
医学肿瘤学临床试验
2019-01-25 10:39圈主:bma 帖子:58
也不是说没用,只是有很多我们不常研究的染色体(如我们研究的主要是常染色体和性染色体以及线粒体,其他的染色体是在染色体间区存在的碱基汇总为一类染色体)会夹杂在我们下载的基因组和基因注释上,将这些不常研究的删除也能减少我们的运算量。你的tophat报错具体报错是怎样的?把整条错误贴过来看下,一般你设置好了参数,把该有的注释和index给tophat是不会报错的,也把你跑的命令发过来给我看下,看下哪里有错。