根据已知的参考基因组和GTF创建适合自己分析的基因组和GTF

  • A+
所属分类:Bioinformatics

###############################################################

######################## Genome build #########################

###############################################################

###删除染色体号后面的内容(包括空格)

###写了一个perl程序,可以切除匹配项之后的信息,得到匹配上的项。

###Perl 名 :cut_space.pl,以下就是程序内的语句:(不同文件需要修改染色体号)

#!/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:

#!/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

#!/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

#!/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

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的方法相同。(略)

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

目前评论:4   其中:访客  4   博主  0

    • avatar zero_jason 1

      删除没用的染色体,怎么判断有用没用的。我运行tophat后的错误是:好象是transcripts有重复。不知道怎么处理呢。求教~~~

        • avatar hygine 1

          @zero_jason 也不是说没用,只是有很多我们不常研究的染色体(如我们研究的主要是常染色体和性染色体以及线粒体,其他的染色体是在染色体间区存在的碱基汇总为一类染色体)会夹杂在我们下载的基因组和基因注释上,将这些不常研究的删除也能减少我们的运算量。
          你的tophat报错具体报错是怎样的?把整条错误贴过来看下,一般你设置好了参数,把该有的注释和index给tophat是不会报错的,也把你跑的命令发过来给我看下,看下哪里有错。

            • avatar zero_jason 1

              @hygine 我的参考基因组名称为,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]

                • avatar baker 1

                  @zero_jason 你是怎么解决的啊,我也遇到这个问题