搜索
楼主: Jimmy

生信编程直播第三题:hg38每条染色体基因,转录本的分布

  [复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-2-5 22:40:32 | 显示全部楼层
huamixinhua 发表于 2017-1-17 18:56
103-python-如梦

[mw_shl_code=python,true]def tj(lj):

你就简单的列了一个代码,你QQ联系我是想说什么呢?我在澳门,所以QQ不是很方便,你加我微信吧,ULWVFJE
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-2-6 19:25:04 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-2-6 19:28 编辑

123-python-qin
计算每条染色体上基因

查看文件数据类型

less -S Homo_sapiens.GRCh38.87.gtf

要求:计算每条染色体上基因数目

就是将GTF文件的第一列chr、第二列中是gene提取出来,将每条染色体上的gene个数求和即可。
[AppleScript] 纯文本查看 复制代码
import re 
import os 
from collections import OrderedDict 
from operator import itemgetter 
os.chdir('/Users/yangqin/Desktop/biotree/python/') 
chr_list = [] #建立数组 
feature_list = {} #建立字典 
with open('Homo_sapiens.GRCh38.87.gtf', 'rt') as f: 
for line in f: 
if line.startswith('#'): #去除注释行 
continue 
line = line.rstrip() #去除右侧换行符 
lst = line.split('\t') #分割数据 
chr_id = lst[0] 
feature = lst[2] 
if chr_id not in chr_list: 
chr_list.append(chr_id) 
feature_list[chr_id] = OrderedDict() #之前没有加这一步,总是出错,参照了别的同学的作业,加上久没有问题了,但是不太明白
feature_list[chr_id][feature] = 0 #在这个字典中,chr_id是key,feature是value 
if feature == "gene": 
feature_list[chr_id][feature] +=1 
for chr_id, feature in feature_list.items(): 
#print (chr_id) 
#print (feature) 
for gene,count in feature.items(): 
print("\t".join([chr_id, str(count)]) + "\n") #输出时,要将数值转换为字符 

结果 显示:
1 5194
2 3971
3 3010
4 2505
5 2868
6 2863
7 2867
X 2359
8 2353
9 2242
11 3235
10 2204
12 2940
13 1304
14 2224
15 2152
16 2511
17 2995
18 1170
20 1386
19 2926
Y 523
22 1318
21 835
MT 37
KI270728.1 8
KI270727.1 8
KI270442.1 2
GL000225.1 1
GL000009.2 1
GL000194.1 2
GL000205.2 2
GL000195.1 3
KI270733.1 4
GL000219.1 1
GL000216.2 1
KI270744.1 1
KI270734.1 4
GL000213.1 2
GL000220.1 4
GL000218.1 1
KI270731.1 2
KI270750.1 1
KI270721.1 4
KI270726.1 2
KI270711.1 1
KI270713.1 4
回复 支持 反对

使用道具 举报

0

主题

13

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-2-7 19:46:56 | 显示全部楼层
本帖最后由 tingting 于 2017-2-7 19:48 编辑

121-R&python-婷

新手上路,还是先完成基础作业再说,高级作业要等把所有作业补完再来做。

第三题比之前两题又有了一些进步,虽然也先看了dongye老师的视频,但是觉得基础作业的题目相对简单,所以尝试着自己写了代码。人类基因组的gtf太大,计算机性能一般,所以就用拟南芥的基因组算了一把,脚本应该大同小异,拟南芥只有5条染色体,加上Pt和Mt,总共7条染色体序列信息。

[Python] 纯文本查看 复制代码
#题目:对人类的基因组注释文件(gtf/gff格式)进行行统计,每条染色体有多少个基因?基因平均有多少个转录本?转录本平均有多个exon?

#思路:
#1.跳过#开头的几行
#2.提取第一列(染色体号),第三列(注释类型)。




for chr_id in ['1','2','3','4','5','Mt','Pt']:
    gene_num = 0
    transcript_num = 0
    exon_num = 0
    print chr_id
    with open('F:test_data/Arabidopsis_thaliana.TAIR10.33.gtf') as f:
        for line in f:
            if line.startswith('#'):
                continue
            line = line.rstrip()
            lst = line.split('\t')
            if lst[0] == chr_id:
                #print lst[0]
                if lst[2] == 'gene':                        
                    #print lst[0],lst[2]
                    gene_num += 1
                elif lst[2] == 'transcript':
                    #print lst[0],lst[2]
                    transcript_num +=1
                elif lst[2] == 'exon':
                    #print lst[0],lst[2]
                    exon_num +=1
                
    print 'chr: %s'%(chr_id)
    print 'gene numbers: %s'%(gene_num)
    print 'transcripts per gene: %.2f' %(transcript_num*1.0/gene_num)
    print 'exons per transcript: %.2f' %(exon_num*1.0/transcript_num)
        
        
        




跑出来的结果是:
1
chr: 1
gene numbers: 8656
transcripts per gene: 1.65
exons per transcript: 5.96
2
chr: 2
gene numbers: 5174
transcripts per gene: 1.65
exons per transcript: 5.65
3
chr: 3
gene numbers: 6435
transcripts per gene: 1.64
exons per transcript: 5.72
4
chr: 4
gene numbers: 4921
transcripts per gene: 1.67
exons per transcript: 5.87
5
chr: 5
gene numbers: 7362
transcripts per gene: 1.64
exons per transcript: 5.91
Mt
chr: Mt
gene numbers: 152
transcripts per gene: 1.00
exons per transcript: 1.12
Pt
chr: Pt
gene numbers: 133
transcripts per gene: 1.00
exons per transcript: 1.18


等下默写完dongye老师的代码再补充到下面。

回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-2-8 05:47:27 | 显示全部楼层
205 perl

目前只重复了视频里的代码,因为太蠢了,抄错了Homo_sapiens.GRCh38.87.gtf中的一个字母,自己debug了一两周都没看见
[AppleScript] 纯文本查看 复制代码
#! usr/bin/perl -w
use strict;

my ( @cells, $cells, @search_feature, $search_feature, @search_cells, @chr, $count, %count, $chr);

open my $row, "Homo_sapiens.GRCh38.87.chr.gtf" or die "Cannot open it";

open OUTPUT,">gene_number.txt";

while (<$row>) {
    last unless $row =~ /\S/;
    @cells = split /\t/, $_;
    if ($cells[2] eq "gene"){
    push @search_feature, $_;
    }
}

foreach $search_feature(@search_feature){
    last unless $search_feature =~ /\S/;
    @search_cells = split /\t/, $search_feature;
    @chr = $search_cells[0];


    foreach $chr(@chr){
    $count{$chr}+=1;
    }
}

foreach $chr(sort keys %count){
    print OUTPUT "$chr\t => \t$count{$chr}\n";
}


结果是这样的
11       =>     3235
12       =>     2940
13       =>     1304
14       =>     2224
15       =>     2152
16       =>     2511
17       =>     2995
18       =>     1170
19       =>     2926
2        =>     3971
20       =>     1386
21       =>     835
22       =>     1318
3        =>     3010
4        =>     2505
5        =>     2868
6        =>     2863
7        =>     2867
8        =>     2353
9        =>     2242
MT       =>     37
X        =>     2359
Y        =>     523

我还会继续修改作业的。蠢死了,耽误好多时间。
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-8 10:24:12 | 显示全部楼层
本帖最后由 chapman 于 2017-2-8 10:28 编辑

[Python] 纯文本查看 复制代码
import re
import sys

args = sys.argv

class Geneome_info():
    def __init__(self):
        self.chr = ""
        self.start = 0
        self.end =0

class Gene(Geneome_info):
    def __init__(self):
        Geneome_info.__init__(self)
        self.orientation = ""
        self.id = ""

class Transcipt(Geneome_info):
    def __init__(self):
        Geneome_info.__init__(self)
        self.id = ""
        self.parent = ""

class Exon(Geneome_info):
    def __init__(self):
        Geneome_info.__init__(self)
        self.parent = ""

        
        
def main(args):
    list_chr = []
    dict_gene = {}
    dict_transcript = {}
    list_exon = []
     
    f = open('Homo_sapiens.GRCh38.87.chr.gtf','r')

    for line in f:
        if line.startswith('#'):
            continue
        lines = line.strip('\n').split('\t')
        Chr = lines[0]
        Type = lines[2]
        start =int(lines[3])
        end = int(lines[4])
        orientation = lines[6]
        attr = lines[8]
        if not re.search(r'protein_coding',attr):
            continue
        if not Chr in list_chr:
            list_chr.append(Chr)

        if Type == "gene":
            gene = Gene()
            Id = re.search(r'gene_id "([^;]+)";?',attr).group(1)
            gene.chr = Chr
            gene.start = start
            gene.end = end
            gene.id = Id
            gene.orientation = orientation
            dict_gene[Id] = gene

        elif Type == "transcript":
            transcript = Transcipt()
            Id = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
            parent = re.search(r'gene_id "([^;]+)";?',attr).group(1)
            if not parent in dict_gene:
                continue
            transcript.chr = Chr
            transcript.start = start
            transcript.end = end
            transcript.id = Id
            transcript.parent = parent
            dict_transcript[Id] = transcript

        elif Type == "exon":
            exon = Exon()
            parent = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
            if not parent in dict_transcript:
                continue
            exon.chr = Chr
            exon.start = start
            exon.end = end
            exon.parent = parent
            list_exon.append(exon)

    
        
    chr_gene(dict_gene)
    gene_len(dict_gene)
    dict_transcript(dict_transcript)
    transcript_exon(list_exon)
    exon_pos(list_exon)

def chr_gene(dict_gene):#染色体上基因数量分布
    count_gene = {}
    for info in dict_gene.values():
        chr = info.chr
        if chr in count_gene:
            count_gene[info.chr] += 1
        else:
            count_gene[info.chr] = 1
    with open('chr_gene.txt','w') as ft_out:
        for che , num in count_gene.items():
            ft_out.write('\t'.join([Chr,str(num)]) + '\n')

            
   
def gene_len(dict_gene):
    with open('gene_len.txt','w') as ft_out:
        for gene_id , info in dict_gene.items():
            Len = info.end - info.start + 1
            ft_out.write('\t'.join([gene_id,str(Len)]) + "\n")
            
def gene_transcript(dict_transcript):
    count_transcript = {}
    for info in list_transcript.values():
        gene_id = info.parent
        if gene_id in count_transcript:
            count_transcript[gene_id] += 1
        else:
            count_transcript[gene_id] = 1
    
    with open('gene_transcript.txt','w') as ft_out:
        for gene_id , num in count_transcript.items():
            Len = info.end - info.start + 1
            ft_out.write('\t'.join([gene_id,str(num)]) + "\n")

def exon_pos(list_exon):
    count_exon = {}
    for exon in list_exon:
        transcript_id =exon.parent
        if transcript_id in count_exon:
            count_exon[transcript_id] += ",%s - %s" % (str(exon.start)),str(exon.end)
        else:
            count_exon[transcript_id] = "%s - %s" % (str(exon.start)),str(exon.end)
    with open('exon_pos.txt','w') as ft_out:
        for transcript_id ,pos in count_exon.items():
            ft_out.write('\t'.join([transcript_id,pos]) + '\t')

#def gene_exon_pos(list_gene , list_transcript,list_exon):
if __name__ =="__main__":
    main(args)

    
    
        

050-Python-晁帅
1、这一讲有点难度,忙活了一个星期才勉强看懂,懂到会写还有一定距离;
2、不熟悉的主要是Python的类和对象不熟悉;
3、如果你也有类似困惑,可以来交流一下。
回复 支持 反对

使用道具 举报

0

主题

12

帖子

165

积分

注册会员

Rank: 2

积分
165
发表于 2017-2-8 20:57:40 | 显示全部楼层
107-perl-Dorom
小白一只,只能跟着视频里的依葫芦画瓢了,因为手头上有大豆的数据,就用大豆的GFF文件做的,只是统计了大豆的每条染色体的基因数,后续的更深层次的问题还有待进一步研究

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
 
my($row,@match,@new_match,@chr,$chr,$count,%count);
open DATA, "Gmax_275_v2.0.gff" or die "Cannot open it";
open OUTPUT,">gene_number.txt";

while($row=<DATA>){
    last unless $row =~ /\S/;
    @match = split /\t/,$row;
    if ($match[2] eq "gene"){
        @new_match = split /\t/,$row;
        @chr = $new_match[0];
            foreach $chr(@chr){
                $count{$chr}+=1;
            }
      }
}
foreach $chr(sort keys %count){
    print OUTPUT "$chr\t=>\t$count{$chr}\n";


运行结果如下:
Chr01   =>      2457
Chr02   =>      3123
Chr03   =>      2649
Chr04   =>      2574
Chr05   =>      2491
Chr06   =>      3258
Chr07   =>      2743
Chr08   =>      3679
Chr09   =>      2865
Chr10   =>      2989
Chr11   =>      2570
Chr12   =>      2425
Chr13   =>      3730
Chr14   =>      2245
Chr15   =>      2774
Chr16   =>      2223
Chr17   =>      2627
Chr18   =>      3023
Chr19   =>      2642
Chr20   =>      2502
回复 支持 反对

使用道具 举报

0

主题

4

帖子

164

积分

注册会员

Rank: 2

积分
164
发表于 2017-2-8 22:49:51 | 显示全部楼层
本帖最后由 DAVIDWANG 于 2017-2-8 22:58 编辑

197R-python davidwang
我侧重真菌基因组的研究,因此对真菌的gff文件进行分析,正则还不是特别懂,就直接拆分选择确定基因id
1 首先是gff文件格式:
##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build MG8
#!genome-build-accession NCBI_Assembly:GCA_000002495.2
##sequence-region CM001231.1 1 7978604
##species https://www.ncbi.nlm.nih.gov/Tax ... wwtax.cgi?id=242507
CM001231.1        Genbank        region        1        7978604        .        +        .        ID=id0;Dbxref=taxon:242507;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA;old-name=Magnaporthe grisea 70-15;strain=70-15
CM001231.1        Genbank        gene        10066        10137        .        -        .        ID=gene0;Name=MGG_20048;gbkey=Gene;gene_biotype=tRNA;locus_tag=MGG_20048;pseudo=true
CM001231.1        Genbank        tRNA        10066        10137        .        -        .        ID=rna0arent=gene0;gbkey=tRNA;product=tRNA-OTHER
CM001231.1        Genbank        exon        10066        10137        .        -        .        ID=id1arent=rna0;gbkey=tRNA;product=tRNA-OTHER
CM001231.1        Genbank        gene        11835        13843        .        -        .        ID=gene1;Name=MGG_01947;end_range=13843,.;gbkey=Gene;gene_biotype=protein_coding;locus_tag=MGG_01947;partial=true;start_range=.,11835
CM001231.1        Genbank        mRNA        11835        13843        .        -        .        ID=rna1arent=gene1;end_range=13843,.;gbkey=mRNA;partial=true;product=cytochrome P450 3A19;start_range=.,11835
CM001231.1        Genbank        exon        13733        13843        .        -        .        ID=id2arent=rna1;end_range=13843,.;gbkey=mRNA;partial=true;product=cytochrome P450 3A19
CM001231.1        Genbank        exon        13399        13468        .        -        .        ID=id3arent=rna1;gbkey=mRNA;partial=true;product=cytochrome P450 3A19
CM001231.1        Genbank        exon        12926        13330        .        -        .        ID=id4arent=rna1;gbkey=mRNA;partial=true;product=cytochrome P450 3A19
CM001231.1        Genbank        exon        12665        12836        .        -        .        ID=id5arent=rna1;gbkey=mRNA;partial=true;product=cytochrome P450 3A19
CM001231.1        Genbank        exon        11835        12558        .        -        .        ID=id6arent=rna1;gbkey=mRNA;partial=true;product=cytochrome P450 3A19;start_range=.,11835
...
其中层级关系:genome<gene<tRNA、mRNA<exon、CDS,每个类型都有自己的id,同时包含蛋白id,也可以通过判断蛋白质的id来确定是否是编码蛋白

代码如下:
[Python] 纯文本查看 复制代码
import sys

args = sys.argv


#根据gff构建类
class GenomeInfo:
    def __init__(self):
        self.chr = ""
        self.start = 0
        self.end = 0

class Gene(GenomeInfo):
    def __init__(self):
        GenomeInfo.__init__(self)
        self.orientation = ""
        self.id = ""

# class Trna(GenomeInfo):
#     def __init__(self):
#         GenomeInfo.__init__(self)
#         self.id = ""
#         self.parent = ""
#
# class Mrna(GenomeInfo):
#     def __init__(self):
#         GenomeInfo.__init__(self)
#         self.id = ""
#         self.parent = ""
#
# class Exon(GenomeInfo):
#     def __init__(self):
#         GenomeInfo.__init__(self)
#         self.id = ""
#         self.parent = ""
#
# class CDS(GenomeInfo):
#     def __init__(self):
#         GenomeInfo.__init__(self)
#         self.id = ""
#         self.parent = ""


def main(args):
    list_chr = []
    list_gene = {}
    with open(args[1]) as fp_gtf:
        for line in fp_gtf:
            if line.startswith("#"):
                continue
            lines = line.strip("\n").split("\t")
            chr = lines[0]
            type = lines[2]
            start = int(lines[3])
            end = int(lines[4])
            orientation = lines[6]
            attr = lines[8]
            if not chr in list_chr:
                list_chr.append(chr)
            if type == "gene":
                ls = attr.split(";")
                id = ls[0].split("=")[1]#利用分割处理得到id
                gene = Gene()
                gene.chr = chr
                gene.start = start
                gene.end = end
                gene.id = id
                gene.orientation = orientation
                list_gene[id] = gene
    chr_gene(list_gene)

def chr_gene(list_gene):
    """
    染色体上基因数量的分布
    :param list_gene:
    :return:
    """
    print("染色体上基因数量的分布:")
    count_gene = {}
    for info in list_gene.values():
        chr = info.chr
        if chr in count_gene:
            count_gene[info.chr] += 1
        else:
            count_gene[info.chr] = 1
    with open ("chr_gene.txt",'w') as fp_out:
        for chr, num in count_gene.items():
            print("\t".join([chr, str(num)])+"\n")
            fp_out.write("\t".join([chr, str(num)])+"\n")


if __name__ == '__main__':
    main(args)



结果如下:
染色体上基因数量的分布:
JH165204.1        1

JH165212.1        1

JH165202.1        1

JH165213.1        1

JH165203.1        1

CM001237.1        1087

JH165184.1        1

CM001231.1        2513

JH165178.1        62

JH165215.1        1

JH165211.1        1

JH165201.1        1

JH165185.1        1

JH165177.1        6

CM001234.1        1830

CM001236.1        1328

JH165175.1        49

JH165218.1        1

CM001235.1        1498

JH165176.1        10

JH165214.1        1

JH165186.1        1

JH165209.1        1

JH165183.1        1

CM001233.1        2149

JH165210.1        1

JH165206.1        1

JH165205.1        1

JH165208.1        1

JH165207.1        1

JH165216.1        1

JH165220.1        1

JH165219.1        1

JH165182.1        1

JH165200.1        1

CM001232.1        2627
看了dongye老师的视频,结合以前对python的了解,有了更深的认识:
1 首先利用class构建了继承基因组类别的几种类,以便容易理解gtf文件的数据结构及继承的信息,同时便于代码重用;
2 利用函数定义数据结构,通过构建的类对数据的结构进行剖析,以便后续的数据个性化分析,逻辑更清晰;
3 练习代码没有捷径,只有反复的敲代码才有感觉。
回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-2-9 17:05:25 | 显示全部楼层
本帖最后由 Aiyawq 于 2017-2-9 17:10 编辑

207-perl-wangQ;
原本我是准备用next if 来进行一行一行筛选的,但是这样做效率太低,于是花了点时间来理解anlan的脚本,终于摸透了二维哈希,这个确实能够节省好多时间呢,目前只能做出这个低级的了解gtf文件的信息,后面进阶了就开始高级任务吧!
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;
open A,$ARGV[0] or die;
#open B,">",$ARGV[1] or die;
my ($a,%hash1,%hash2);
print "chr\tgene\tprotein_coding\ttranscript\texon\n";
while(<A>){
       chomp;
       next if /^#/;
       my @a = split (/\t/, $_);
       $hash1{$a[0]}{$a[2]}++;
       if($a[2] eq "gene"){
               if($a[8] =~ /.*?(protein_coding)/){
                       $hash1{$a[0]}{$1}++;
               }
       }
}
close A;
foreach(sort keys %hash1){
        print"$_\t$hash1{$_}->{gene}\t$hash1{$_}->{protein_coding}\t$hash1{$_}->{transcript}\t$hash1{$_}->{exon}\n";
}
#close B;
结果如下:

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

2

帖子

16

积分

新手上路

Rank: 1

积分
16
发表于 2017-2-9 17:40:47 | 显示全部楼层
[AppleScript] 纯文本查看 复制代码
source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(GenomicFeatures) #根据Google搜索得出
setwd("E:/BaiduNetdiskDownload") #下载Homo_sapiens.GRCh38.87.chr.gtf.gz的存放地点
GRCh38_txdb<-makeTxDbFromGFF("Homo_sapiens.GRCh38.87.chr.gtf.gz",organism = "Homo sapiens")
columns(GRCh38_txdb) #染色体例的信息
seqlevels(GRCh38_txdb) #多少染色体

#查看各条染色体的基因分布
grch38_gene<-genes(GRCh38_txdb,columns="gene_id")
grch38_gene_df<-as.data.frame(grch38_gene)
head(grch38_gene_df) 
nrow(grch38_gene_df) #总基因的个数
#[1] 57992
library(ggplot2)
#geom_bar是对数目的作图
ggplot(grch38_gene_df,aes(x=factor(seqnames)))+geom_col()

#同理查看各转录本
grch38_tx<-transcripts(GRCh38_txdb,columns="tx_name")
grch38_tx_df<-as.data.frame(grch38_tx)
nrow(grch38_tx_df)
#[1] 197935
ggplot(grch38_tx_df,aes(x=factor(seqnames)))+geom_bar()#得到转录本按染色体的分布
 
#在查看各个基因所含转录本的平均水平
GRCh38_tx_by_genne<-transcriptsBy(GRCh38_txdb,by=c("gene"))
GRCh38_tx_by_genne<- as.data.frame(GRCh38_tx_by_genne)
nrow(GRCh38_tx_by_genne)
#[1] 197935
library(dplyr)
GRCh38_tx_by_genne<-GRCh38_tx_by_genne%>%group_by(group_name)%>%mutate(count=n())
head(GRCh38_tx_by_genne) #基因出现的次数,转录本
sort(as.numeric(unique(GRCh38_tx_by_genne$count)))
summary(GRCh38_tx_by_genne$count) 对COUNT取平均数as.numeric
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1.00    3.00    9.00   11.04   15.00  170.00
df<-unique(GRCh38_tx_by_genne[,c("group_name","count")]) 具体分布情况,基因中所含转录本数
ggplot(df,aes(x=factor(count)))+geom_bar()

#考察转录本所含外显子数目问题
GRCh38_tx_by_exon<-transcriptsBy(GRCh38_txdb,by="exon")
GRCh38_tx_by_exon<- as.data.frame(GRCh38_tx_by_genne)
nrow(as.data.frame(GRCh38_tx_by_exon))
#[1] 1181908
GRCh38_tx_by_exon<-as.data.frame(GRCh38_tx_by_exon)
GRCh38_tx_by_exon<-GRCh38_tx_by_exon%>%group_by(tx_name)%>%mutate(count=n())
df_exon<-unique(GRCh38_tx_by_exon[,c("tx_name","count")])
summary(df_exon)
group_name            count        
Length:57992       Min.   :  1.000  
Class :character   1st Qu.:  1.000  
Mode  :character   Median :  1.000  
                    Mean   :  3.413  
                    3rd Qu.:  4.000  
                    Max.   :170.000  
ggplot(df_exon,aes(x=factor(count)))+geom_bar()
回复 支持 反对

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-9 21:30:30 | 显示全部楼层
本帖最后由 zckoo007 于 2017-2-9 21:45 编辑

002+Python+R

1. 下载pycharm  摸索了三天,了解一些基本功能,project创建(pycharm 相当强大!)
2. 下载Git bash  一天,应为不会Linux,以后慢慢摸索
3. python  函数的定义 两天,
4. python class 一天
5. 读懂代码 两天
6. 感谢 dongye 哥!
[Python] 纯文本查看 复制代码
import sys
import re

args = sys.argv
'''sys.argv 是命令行参数,是一个字符串列表,0代表其路径'''


class Genome_info:      #这是一个基类,所有类型都通用的
    def __init__(self):
        self.chr = ""   #染色体号
        self.start = 0
        self.end = 0


class Gene(Genome_info):   #这个函数继承了基类
    def __init__(self):
        Genome_info.__init__(self)
        self.orientation = ""
        self.id = ""


class Transcript(Genome_info):
    def __init__(self):
        Genome_info.__init__(self)
        self.id = ""
        self.parent = ""


class Exon(Genome_info):
    def __init__(self):
        Genome_info.__init__(self)
        self.parent = ""


def main(args):
    """
    一个输入参数:
    第一个参数为物种gtf文件

    :return:
    """

    list_chr = []
    list_gene = {}   #因为有 id 号 所以用dir
    list_transcript = {}
    list_exon = []
    # l_n = 0
    with open(args[1]) as fp_gtf:        #打开文件遍历GTF
        for line in fp_gtf:
            if line.startswith("#"):      #号开头的注释过滤掉,不统计这个
                continue
            # print ("in %s" % l_n)
            # l_n += 1
            lines = line.strip("\n").split("\t")
            chr = lines[0]
            type = lines[2]
            start = int(lines[3])
            end = int(lines[4])
            orientation = lines[6]
            attr = lines[8]
            if not re.search(r'protein_coding', attr):   #取到的是蛋白的编码,如果没有的话就跳过,不做统计了
                continue

            if not chr in list_chr:    #把染色体的类型添加到列表里
                list_chr.append(chr)

            if type == "gene":
                gene = Gene()          #初始化一个基因对象
                id = re.search(r'gene_id "([^;]+)";?', attr).group(1)  # 0 返回所有列表,1取第一个
                gene.chr = chr
                gene.start = start
                gene.end = end
                gene.id = id
                gene.orientation = orientation
                list_gene[id] = gene
                # print(id)
            elif type == "transcript":
                transcript = Transcript()
                id = re.search(r'transcript_id "([^;]+)";?', attr).group(1)
                parent = re.search(r'gene_id "([^;]+)";?', attr).group(1)
                if not parent in list_gene:
                    continue
                transcript.chr = chr
                transcript.start = start
                transcript.end = end
                transcript.id = id
                transcript.parent = parent
                list_transcript[id] = transcript

            elif type == "exon":
                exon = Exon()
                parent = re.search(r'transcript_id "([^;]+)";?', attr).group(1)
                if not parent in list_transcript:
                    continue
                exon.chr = chr
                exon.start = start
                exon.end = end
                exon.parent = parent
                list_exon.append(exon)

    chr_gene(list_gene)
    gene_len(list_gene)
    gene_transcript(list_transcript)
    transcript_exon(list_exon)
    exon_pos(list_exon)


def chr_gene(list_gene):
    """
    染色体上基因数量分布

    :param list_gene:
    :return:
    """

    print("染色体上基因数量分布")
    count_gene = {}     #这是一个计数器
    for info in list_gene.values():
        chr = info.chr
        if chr in count_gene:
            count_gene[info.chr] += 1
        else:
            count_gene[info.chr] = 1
    with open("chr_gene.txt", 'w') as fp_out:
        for chr, num in count_gene.items():
            print("\t".join([chr, str(num)]) + "\n")
            fp_out.write("\t".join([chr, str(num)]) + "\n")


def gene_len(list_gene):
    """
    基因长度分布情况

    :param list_gene:
    :return:
    """

    print("基因长度分布情况")
    with open("gene_len.txt", 'w') as fp_out:
        for gene_id, info in list_gene.items():
            len = info.end - info.start + 1
            fp_out.write("\t".join([gene_id, str(len)]) + "\n")
            print("\t".join([gene_id, str(len)]) + "\n")


def gene_transcript(list_transcript):
    """
    基因的转录本数量分布

    :param list_transcript:
    :return:
    """

    print("基因的转录本数量分布")
    count_transcript = {}
    for info in list_transcript.values():
        gene_id = info.parent
        if gene_id in count_transcript:
            count_transcript[gene_id] += 1
        else:
            count_transcript[gene_id] = 1
    with open("gene_transcript.txt", 'w') as fp_out:
        for gene_id, num in count_transcript.items():
            print("\t".join([gene_id, str(num)]) + "\n")
            fp_out.write("\t".join([gene_id, str(num)]) + "\n")


def transcript_exon(list_exon):
    """
    转录本的外显子数量统计

    :param list_exon:
    :return:
    """

    print("转录本的外显子数量统计")
    count_exon = {}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id in count_exon:
            count_exon[transcript_id] += 1
        else:
            count_exon[transcript_id] = 1
    with open("transcript_exon.txt", 'w') as fp_out:
        for transcript_id, num in count_exon.items():
            print("\t".join([transcript_id, str(num)]) + "\n")
            fp_out.write("\t".join([transcript_id, str(num)]) + "\n")


def exon_pos(list_exon):
    """
    外显子坐标统计

    :param list_exon:
    :return:
    """

    print("外显子坐标统计")
    count_exon = {}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id in count_exon:
            count_exon[transcript_id] += ",%s-%s" % (str(exon.start), str(exon.end))
        else:
            count_exon[transcript_id] = "%s-%s" % (str(exon.start), str(exon.end))
    with open("exon_pos.txt", 'w') as fp_out:
        for transcript_id, pos in count_exon.items():
            print("\t".join([transcript_id, pos]) + "\n")
            fp_out.write("\t".join([transcript_id, pos]) + "\n")


def gene_exon_pos(list_gene, list_transcript, list_exon):
    """
    根据exon的parent将所有exon对应到transcript
    根据transcript的parent将所有transcript对应到gene
    根据gene按chr分组得到chromosome列表

    从chromosome中输出某个指定基因的所有外显子坐标信息并画图
    生信编程直播第五题

    :param list_gene:
    :param list_transcript:
    :param list_exon:
    :return:
    """
    pass


if __name__ == "__main__":
    main(args)



染色体上基因数量分布
1        5

基因长度分布情况
ENSG00000186092        918

ENSG00000279928        1766

ENSG00000279457        15400

ENSG00000278566        939

ENSG00000273547        939

基因的转录本数量分布
ENSG00000186092        1

ENSG00000279928        1

ENSG00000279457        3

ENSG00000278566        1

ENSG00000273547        1

转录本的外显子数量统计
ENST00000335137        1

ENST00000624431        3

ENST00000623834        10

ENST00000623083        11

ENST00000624735        15

ENST00000426406        1

ENST00000332831        1

外显子坐标统计
ENST00000335137        69091-70008

ENST00000624431        182393-182746,183114-183240,183922-184158

ENST00000623834        195259-195411,188791-188892,188489-188584,188439-188486,188130-188266,188022-188028,187755-187886,187376-187577,185491-185559,184923-185350

ENST00000623083        195263-195411,188791-188902,188489-188584,188439-188486,188126-188266,187755-187886,187376-187577,187129-187287,186317-186469,185491-185559,184925-185350

ENST00000624735        200050-200322,195259-195416,188791-188889,188489-188584,188439-188486,188130-188266,188102-188105,187755-187886,187380-187577,187270-187287,187129-187267,186317-186469,185529-185559,184977-185049,184927-184971

ENST00000426406        450740-451678

ENST00000332831        685716-686654


Process finished with exit code 0
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 15:39 , Processed in 0.061659 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.