搜索
查看: 622|回复: 19

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

[复制链接]

276

主题

477

帖子

1729

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1729
发表于 2017-1-5 08:36:03 | 显示全部楼层 |阅读模式
这一讲主要针对,gtf注释文件的探究!
我写过一个R完成的版本:R的bioconductor包TxDb.Hsapiens.UCSC.hg19.knownGene详解 :http://www.bio-info-trainee.com/831.html
人类的基因组注释文件,一般是gtf/gff格式的。参考:http://www.bio-info-trainee.com/782.html
用perl和python,R,shell均可以完成!基础作业,就是对这个文件ftp://ftp.ensembl.org/pub/releas ... RCh38.87.chr.gtf.gz  进行统计,普通人只需要关心第1,3列就好了。
每条染色体基因个数的分布?所有基因平均有多少个转录本?所有转录本平均有多个exon和intron?或者其它一系列你想统计的,你想了解的!
或者其它一系列你想统计的,你想了解的!
或者其它一系列你想统计的,你想了解的!


比如我们先看每条染色体的基因个数:
[Shell] 纯文本查看 复制代码
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{print if $F[2] eq "gene" }' |cut -f 1 |sort |uniq -c

   5194 1
   2204 10
   3235 11
   2940 12
   1304 13
   2224 14
   2152 15
   2511 16
   2995 17
   1170 18
   2926 19
   3971 2
   1386 20
    835 21
   1318 22
   3010 3
   2505 4
   2868 5
   2863 6
   2867 7
   2353 8
   2242 9
     37 MT
   2359 X
    523 Y
但是如果只考虑protein coding gene的话:
2052 1
    732 10
   1276 11
   1034 12
    327 13
    623 14
    597 15
    866 16
   1197 17
    270 18
   1470 19
   1255 2
    544 20
    233 21
    438 22
   1076 3
    751 4
    876 5
   1045 6
    906 7
    676 8
    781 9
     13 MT
    841 X
     53 Y
Y染色体的基因数量明显减少了,可以看到,很多基因都是其它类型的!
那么基因分类是怎么样的呢?
[Shell] 纯文本查看 复制代码
 zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}'  |sort |uniq -c

   30 3prime_overlapping_ncRNA
   5526 antisense
      4 bidirectional_promoter_lncRNA
     14 IG_C_gene
      9 IG_C_pseudogene
     37 IG_D_gene
     18 IG_J_gene
      3 IG_J_pseudogene
      1 IG_pseudogene
    145 IG_V_gene
    193 IG_V_pseudogene
   7533 lincRNA
      1 macro_lncRNA
   1567 miRNA
   2212 misc_RNA
      2 Mt_rRNA
     22 Mt_tRNA
      3 non_coding
     51 polymorphic_pseudogene
  10268 processed_pseudogene
    515 processed_transcript
  19932 protein_coding
     21 pseudogene
      8 ribozyme
    543 rRNA
     49 scaRNA
      1 scRNA
    908 sense_intronic
    187 sense_overlapping
    944 snoRNA
   1900 snRNA
      5 sRNA
   1048 TEC
    450 transcribed_processed_pseudogene
     60 transcribed_unitary_pseudogene
    735 transcribed_unprocessed_pseudogene
      6 TR_C_gene
      4 TR_D_gene
     79 TR_J_gene
      4 TR_J_pseudogene
    108 TR_V_gene
     30 TR_V_pseudogene
    154 unitary_pseudogene
   2661 unprocessed_pseudogene
      1 vaultRNA

我希望你们可以发挥创造力,把这个gtf文件玩透彻,彻底了解基因结构。
我希望你们可以发挥创造力,把这个gtf文件玩透彻,彻底了解基因结构。
我希望你们可以发挥创造力,把这个gtf文件玩透彻,彻底了解基因结构。

我只是随便举几个可以统计的例子,要学会生信编程,还是要靠你们自己的。

高级作业,下载human,rat,mouse,dog,cat,chicken,等十几种物种的gtf文件:http://asia.ensembl.org/info/data/ftp/index.html 然后写好染色体基因,转录本的分布统计函数,转录本外显子个数统计分布函数,对它们批量处理。

所有基因平均有多少个转录本?所有转录本平均有多个exon和intron?
如果要比较多个数据库呢?gencode?UCSC?NCBI?
如果把基因分成多个类型呢?protein coding gene,pseudogene,lncRNA还有miRNA的基因?它们的特征又有哪些变化呢?
我希望看到一个类似下面的图:






本帖子中包含更多资源

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

x



上一篇:生信编程直播第二题-hg19基因组序列的一些探究
下一篇:有关Drop-Seq技术的一篇报道
回复

使用道具 举报

0

主题

10

帖子

83

积分

注册会员

Rank: 2

积分
83
发表于 2017-1-15 21:58:56 | 显示全部楼层
本帖最后由 x2yline 于 2017-1-23 17:09 编辑

077-R-python-shell-x2yline

主要文字来源:于浩,张晶,张英主编.生物信息学 实验指导 第1版[M].长春:吉林大学出版社.2009.

Ensemble( ensembl.org网站是常用真核生物参考基因组来源之一 )能够对人类基因自动进行注释,包括人类,小鼠,斑马鱼,猪和大鼠等,也包括来自HAVANA的人工注释信息。
Ensembl是一项生物信息学研究计划,旨在开发种能够对真核生物基因组进行自动注释(automatic annotation)并加以维护的软件系统。该计划由英国Sanger研究所Wellcome基金会及欧洲分子生物学实验室所属分部欧洲生物信息学研究所共同协作运营。

Ensembl与NCBI的NCBI Map Viewer和UCSC是最为常用基因组检索数据库。

Ensembl 与NCBI Map Viewer和UCSC最大区别表现在以下5点:
a.Ensembl的基因数据集是依据mRNA和蛋内序列的数据信息白动注释的。数据来源为新的基因组数据,UniProt/SwissProt和UniProt/TrEMBL的蛋白序列,NCBI的RefSeq里的DNA和蛋白序列和EMBL的cDNA序列。
b.Ensembl是一个开源(Perl API )的全自动的基因注释软件系统,很多网站都采用Ensembl这套软件系统。
c.Ensembl拥存其特有的BioMart功能。BioMart可以依据设定的要求对基 因组进行条件性检索,检索的结果吋以以图表的形式给出。
d.与其它数据库相整合,比如DAS。
e.基因组间的比较分析。

基因注释机构
目前从事基因注释的机构组织有很多,这里列出的只是较为常用的几个。
1. Ensembl:目的是做出最好的基因注释集。
2.Havana (VEGA):是桑格中心的一个基因注释组织,它的目标和Eiisembl—致,因此,结合得也最紧密。
3. HGNC -给出人类基因唯一的名字和符号。
4. UniProt 主要集中于蛋白质的信息注释。

Ensembl的通用基因注释有两种,一是Ensembl GeneBuild,它是自动化注释,速度快,实时更新,在不同物种上均适用;另一种是Wellcome基金会的 Havana (VEGA)小组的注释,它是手工注释,速度慢,但是准确,它依据的都是已经验证过的mRNA和蛋白序列来注释,比较费时。因此Ensembl基因组数据库 中,会有两种注释。

Havana (VEGA)小组的注释常有以下几种类型
详细信息:http://vega.sanger.ac.uk/info/about/gene_and_transcript_types.html
Protein coding: 包括开放阅读框 (ORF).
Processed transcript:没有开放阅读框(ORF)
Pseudogene:假基因,是指脱氧核糖核酸(DNA)的碱基序列中,一段与其他生物体内已知的基因序列非常相似的片段。但是这个片段由于移码突变或者无义突变破坏了ORF,无法发挥原有的基因功能,也就是无法制造出蛋白质
IG gene:免疫球蛋白家族基因
TR Gene:T细胞受体基因
TEC (To be Experimentally Confirmed)

人类和小鼠基因组的GTF文件与GENCODE计划发布的gene set文件相同。
The GENCODE project 的目标为对人类和小鼠基因组提供高质量的注释信息和实验确证。
The GENCODE gene sets被其他项目作为参考而广泛使用(如 1000 Genomes).
详细内容:https://www.gencodegenes.org/about.html


带有abinitio扩展名的文件为用Genescan和abinitio基因预测工具生成的
预测基因的注释文件


------------------
GTF文件输出示例
------------------

#!genome-build GRCh38
11      ensembl_havana  gene    5422111 5423206 .       +       .       gene_id "ENSG00000167360"; gene_version "4"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
11      ensembl_havana  transcript      5422111 5423206 .       +       .       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381";
11      ensembl_havana  exon    5422111 5423206 .       +       .       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; exon_number "1"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381"; exon_id "ENSE00001276439"; exon_version "4";
11      ensembl_havana  CDS     5422201 5423151 .       +       0       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; exon_number "1"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381"; protein_id "ENSP00000300778"; protein_version "4";
11      ensembl_havana  start_codon     5422201 5422203 .       +       0       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; exon_number "1"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381";
11      ensembl_havana  stop_codon      5423152 5423154 .       +       0       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; exon_number "1"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381";
11      ensembl_havana  UTR     5422111 5422200 .       +       .       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381";
11      ensembl_havana  UTR     5423155 5423206 .       +       .       gene_id "ENSG00000167360"; gene_version "4"; transcript_id "ENST00000300778"; transcript_version "4"; gene_name "OR51Q1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "OR51Q1-001"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS31381";

1.python实现:
[AppleScript] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

import time
import re
from collections import OrderedDict
start = time.clock()
 
def count_gtf_gene(file_path, buffer_size=1024*1024):
    print(file_path.split('.')[-1])
    gtf_count = OrderedDict()
    with open(file_path, 'r') as f:
        for i in f:
            if not i.startswith('#'):
                i_list = i.split('\t')
                CHR = i_list[0]
                feature = i_list[2]
                if CHR not in gtf_count.keys():
                    gtf_count[CHR] = OrderedDict()
                    gtf_count[CHR]['all_gene'] = 0
                if feature == 'gene':
                    gtf_count[CHR]['all_gene'] += 1
                    TYPE = re.search('gene_biotype "(.*?)";',i_list[-1]).group(1)
                    try:
                        gtf_count[CHR][TYPE] += 1
                    except:
                        gtf_count[CHR][TYPE] = 1
    return gtf_count
 
def write_to_csv(gtf_gene_dict,file_save):   
    c = []
    for a,b in gtf_gene_dict.items():
        for i in b.keys():
            c.append(i)
    c = set(c)
    file_raw = 'Chromsome'
    file_raw += ',' + 'all_gene' + ',' + 'protein_coding'
    for i in c:
        if i not in ['all_gene', 'protein_coding']:
            file_raw += ',' + i
    for chr_id in list(range(1,23))+['X','Y','MT']:
        chr_id = str(chr_id)
        file_raw += '\nchr_'+ chr_id
        file_raw += ',' + str(gtf_gene_dict[chr_id]['all_gene']) +',' + str(gtf_gene_dict[chr_id]['protein_coding'])
        for i in c:
            if i not in ['all_gene', 'protein_coding']:
                i = str(i)
                try:
                    file_raw += ',' + str(gtf_gene_dict[chr_id])
                except:
                    file_raw += ',' + str(0)
    file_raw_list = file_raw.split('\n')[1:]
    i_list = file_raw_list[1].split(',')[1:]
    SUM = list(range(len(i_list)))
    for j in range(len(i_list)):
        SUM[j] = 0
    for i in file_raw_list:
        i_list = i.split(',')[1:]
        for j in range(len(i_list)):
            SUM[j] += int(i_list[j])
    file_raw += '\nSUM'
    for i in SUM:
        file_raw += ',' + str(i)
    with open(file_save,'w') as f:
        f.write(file_raw)
    return re.sub(',', '\t', file_raw)
 
def count_transcript_and_exon(file_path):
    print(file_path)
    transcript = []
    exon = []
    with open(file_path, 'r') as f:
        for i in f:
            if not i.startswith('#'):
                i_list = i.split('\t')
                feature = i_list[2]
                if feature == 'gene':
                    try:
                        transcript.append(transcript_num)
                    except:
                        transcript_num = 0
                    transcript_num = 0
                elif feature == 'transcript':
                    try:
                        exon.append(exon_num)
                    except:
                        exon_num = 0
                    transcript_num += 1
                    exon_num = 0
                elif feature == 'exon':
                    exon_num += 1
        try:
            transcript.append(transcipt_num)
        except:
            pass
        try:
            exon.apend(exon_num)
        except:
            pass
                     
    return (transcript, exon)
 
     
file_path_human = r'F:\genome\Homo_sapiens.GRCh38.87.chr.gtf\Homo_sapiens.GRCh38.87.chr.gtf'
#gtf_gene_dict1 = count_gtf_gene(file_path1)
#file_to_write1 = write_to_csv(gtf_gene_dict, file_save = 'gtf_analysis.csv')
file_path_cat = r'F:\genome\Felis_catus.Felis_catus_6.2.87.chr.gtf\Felis_catus.Felis_catus_6.2.87.chr.gtf'
#print('Used %.4f s'%(time.clock()-start))   
file_path_gorilla = r'F:\genome\Gorilla_gorilla.gorGor3.1.87.chr.gtf\Gorilla_gorilla.gorGor3.1.87.chr.gtf'
file_path_mus = r'F:\genome\Mus_musculus.GRCm38.87.chr.gtf\Mus_musculus.GRCm38.87.chr.gtf'
file_path_zeb = r'F:\genome\Danio_rerio.GRCz10.87.chr.gtf\Danio_rerio.GRCz10.87.chr.gtf'
file_path_chicken = r'F:\genome\Gallus_gallus.Gallus_gallus-5.0.87.chr.gtf\Gallus_gallus.Gallus_gallus-5.0.87.chr.gtf'
file_path_elephant = r'F:\genome\Caenorhabditis_elegans.WBcel235.87.gtf\Caenorhabditis_elegans.WBcel235.87.gtf'
file_path_fluitfly = r'F:\genome\Drosophila_melanogaster.BDGP6.87.chr.gtf\Drosophila_melanogaster.BDGP6.87.chr.gtf'
file_path_pig = r'F:\genome\Sus_scrofa.Sscrofa10.2.87.chr.gtf\Sus_scrofa.Sscrofa10.2.87.chr.gtf'
file_path_dog = r'F:\genome\Canis_familiaris.CanFam3.1.87.chr.gtf\Canis_familiaris.CanFam3.1.87.chr.gtf'
file_path_horse = r'F:\genome\Equus_caballus.EquCab2.87.chr.gtf\Equus_caballus.EquCab2.87.chr.gtf'

transcript_human, exon_human = count_transcript_and_exon(file_path_human)
transcript_cat, exon_cat = count_transcript_and_exon(file_path_cat)
transcript_zeb, exon_zeb = count_transcript_and_exon(file_path_zeb)
transcript_gorilla, exon_gorilla = count_transcript_and_exon(file_path_gorilla)
transcript_mus, exon_mus = count_transcript_and_exon(file_path_mus)
transcript_pig, exon_pig = count_transcript_and_exon(file_path_pig)
transcript_fluitfly, exon_fluitfly = count_transcript_and_exon(file_path_fluitfly)
transcript_dog, exon_dog = count_transcript_and_exon(file_path_dog)
transcript_chicken, exon_chicken = count_transcript_and_exon(file_path_chicken)
transcript_horse, exon_horse = count_transcript_and_exon(file_path_horse)
transcript_elephant, exon_elephant = count_transcript_and_exon(file_path_elephant)
exon_hub = [exon_human, exon_gorilla, exon_elephant, exon_horse, exon_pig, exon_dog, exon_cat,  exon_mus, exon_zeb, exon_fluitfly]

def list_to_freq(exon):
    exon_range = [0]
    exon_freq = [0]
    for i in set(exon):
        exon_range.append(i)
        exon_freq.append(exon.count(i)/len(exon)*100)
    return (exon_range, exon_freq)
 
import numpy as np
import matplotlib
from pylab import *
import matplotlib.pyplot as plt
myfont = matplotlib.font_manager.FontProperties(fname='C:/Windows/Fonts/msyh.ttf') 
mpl.rcParams['axes.unicode_minus'] = False

plt.figure(1)
plt.style.use('seaborn-white')
plt.xlim(0,15)
plt.ylim(0, 50)
plt.ylabel('百分比(%)',fontproperties=myfont)
plt.xlabel(u'外显子数目',fontproperties=myfont)
plt.title(u'不同物种基因外显子数目频率统计图', fontproperties=myfont)
i = -1
colors = ['#00FF00','#006400','blue','#9400D3','#800000','#191970','#CD853F','#FF0000','#FFFF00','#FFA500']
for exon_item in exon_hub:
    i += 1
    labels = ['Human','gorilla','elephant', 'horse', 'pig', 'dog', 'cat',  'mus', 'zebrafish', 'fluitfly' ]
    #colors = ['r', 'g', 'k', 'b','r', 'g', 'k', 'b', 'r', 'b']
    exon_range, exon_freq = list_to_freq(exon_item)
    plt.plot(exon_range, exon_freq, color=colors, alpha=.61, label=labels)
    #plt.xticks(exon_range[0:15])
    plt.plot((median(exon_item)+i/100,median(exon_item)+i/100), (0,50), colors,ls = '--', alpha=1)
    #plt.text(median(exon_item)-0.5+i/10,2,u'中位数', ha='left', color= colors, fontproperties=myfont, alpha=0.8)
 
plt.legend(loc='upper right')
plt.savefig('F:\\1.png', dpi=600)



plt.figure(2)
plt.style.use('ggplot')
plt.boxplot(exon_hub)
plt.ylim(0, 35)
plt.yticks(range(1,36,2))
plt.xticks(range(1,11),labels, fontsize=9)
plt.title('exons of different species')
plt.ylabel('exon number')
plt.savefig('F:\\2.png', dpi=600)
plt.show()











贴上human gtf文件跑出来的部分结果
Chromsome        all_gene        protein_coding        polymorphic_pseudogene        processed_pseudogene        IG_C_pseudogene        scRNA        unitary_pseudogene        TR_J_gene        sense_overlapping        IG_D_gene        misc_RNA
chr_1        5194        2052        6        897        0        0        26        0        17        0        192
chr_2        3971        1255        2        751        0        1        7        0        10        0        176
chr_3        3010        1076        2        616        0        0        9        0        7        0        134
chr_4        2505        751        1        561        0        0        2        0        9        0        104
chr_5        2868        876        0        625        0        0        4        0        10        0        119
chr_6        2863        1045        3        639        0        0        10        0        7        0        105
chr_7        2867        906        2        577        0        0        13        19        12        0        143
chr_8        2353        676        1        468        0        0        4        0        8        0        82
chr_9        2242        781        2        461        1        0        5        0        9        0        96
chr_10        2204        732        1        422        0        0        5        0        9        0        89
chr_11        3235        1276        22        473        0        0        20        0        17        0        97
chr_12        2940        1034        0        490        0        0        6        0        5        0        115
chr_13        1304        327        0        295        0        0        1        0        0        0        75
chr_14        2224        623        3        324        2        0        2        60        11        27        79
chr_15        2152        597        0        287        0        0        1        0        8        10        93
chr_16        2511        866        1        252        0        0        12        0        12        0        51
chr_17        2995        1197        0        310        0        0        6        0        3        0        99
chr_18        1170        270        0        207        1        0        3        0        1        0        41
chr_19        2926        1470        2        263        0        0        7        0        6        0        61
chr_20        1386        544        0        198        0        0        5        0        3        0        68
chr_21        835        233        0        143        0        0        0        0        8        0        24
chr_22        1318        438        2        156        5        0        3        0        9        0        62
chr_X        2359        841        1        723        0        0        3        0        6        0        100
chr_Y        523        53        0        130        0        0        0        0        0        0        7
chr_MT        37        13        0        0        0        0        0        0        0        0        0
SUM        57992        19932        51        10268        9        1        154        79        187        37        2212

最近正好看到有关染色体与基因关系的一页课件,以下步骤只能说是对课件中的结论进行验证
'探索代码如下:
[Python] 纯文本查看 复制代码
import time
import re
from collections import OrderedDict
start = time.clock()
  
def count_gtf_gene(file_path, buffer_size=1024*1024):
    print(file_path.split('.')[-1])
    gtf_count = OrderedDict()
    with open(file_path, 'r') as f:
        for i in f:
            if not i.startswith('#'):
                i_list = i.split('\t')
                CHR = i_list[0]
                feature = i_list[2]
                if CHR not in gtf_count.keys():
                    gtf_count[CHR] = OrderedDict()
                    gtf_count[CHR]['all_gene'] = 0
                if feature == 'gene':
                    gtf_count[CHR]['all_gene'] += 1
                    TYPE = re.search('gene_biotype "(.*?)";',i_list[-1]).group(1)
                    try:
                        gtf_count[CHR][TYPE] += 1
                    except:
                        gtf_count[CHR][TYPE] = 1
    return gtf_count
      
file_path_human = r'F:\genome\Homo_sapiens.GRCh38.87.chr.gtf\Homo_sapiens.GRCh38.87.chr.gtf'

hg38_gene = count_gtf_gene(file_path_human)
chr_id_list = list(range(1,23)) + ['X','Y','MT']
gene_num_list = []
for i in chr_id_list:
    gene_num_list.append(hg38_gene[str(i)]['all_gene'])

import matplotlib.pyplot as plt
plt.figure(1)
plt.style.use('seaborn-white')
plt.bar(left = range(25), height = gene_num_list, color='k')
for i in range(len(gene_num_list)):
    plt.text( x= i , y=gene_num_list+35,s=str(round(gene_num_list)), fontsize = 5)
plt.title('Gene distribution of hg38 genome')
plt.ylabel('Gene number of chromosome')
pos = []
for i in range(len(chr_id_list)):
    pos.append(i + 0.35)
plt.xticks(pos, list(range(1,23)) + ['X','Y','MT'], fontsize=8)
plt.xlim(-0.2, )
plt.savefig('F:\gene_distribution.png',dpi=600)
plt.show()




plt.figure(2)
plt.style.use('seaborn-white')
import os
from collections import OrderedDict

def count_fasta_atcgn(file_path, buffer_size=1024*1024):
    bases = ['N', 'A', 'T', 'C', 'G']
    ATCG_analysis = OrderedDict()
    with open(file_path, 'r') as f:
        line1 = f.readline().upper()
        chr_i = re.split('\s', line1)[0][1:]
        print(chr_i)
        ATCG_analysis[chr_i] = OrderedDict()
        for base in bases:
            ATCG_analysis[chr_i][base] = 0
        while True:
            chunk = f.read(buffer_size).upper()
            if '>' in chunk:
                chromsome = re.split('>',chunk)
                if chromsome[0]:
                    for base in bases:
                        ATCG_analysis[chr_i][base] += chromsome[0].count(base)
                for i in chromsome[1:]:
                    if i:
                        chr_i = re.split('\s', i[0:i.index('\n')])[0]
                        print(chr_i)
                        strings_i = i[i.index('\n'):]
                        ATCG_analysis[chr_i] = OrderedDict()
                        for base in bases:
                            ATCG_analysis[chr_i][base] = strings_i.count(base)
            else:
                for base in bases:
                    ATCG_analysis[chr_i][base] += chunk.count(base)
            if not chunk:
                break
    return ATCG_analysis
 
file_path = r'F:\genome\hg38.chromFa\chroms\hg38.fa'

ATCG_analysis = count_fasta_atcgn(file_path, buffer_size=1024*1024)
cg_list = []
chr_id_list = list(range(1,23)) + ['X','Y','M']
for i in chr_id_list:
    cg_list.append((ATCG_analysis['CHR'+str(i)]['G']+ATCG_analysis['CHR'+str(i)]['C'])/(ATCG_analysis['CHR'+str(i)]['A']+ATCG_analysis['CHR'+str(i)]['T']+ATCG_analysis['CHR'+str(i)]['C']+ATCG_analysis['CHR'+str(i)]['G'])*100)

plt.bar(left = range(25), height = cg_list, color='k')
for i in range(len(cg_list)):
    plt.text( x=i -.1, y=cg_list+.35,s=str(round(cg_list)))
plt.title('GC content for hg38 genome')
plt.ylabel('GC content (%)')
pos = []
for i in range(len(chr_id_list)):
    pos.append(i + 0.35)
plt.xticks(pos, list(range(1,23)) + ['X','Y','MT'], fontsize=8)
plt.xlim(-0.2, )
plt.ylim(0, 100)
plt.savefig('F:\hg38_gc.png',dpi=600)
plt.show()


结果(图)为:





课件如下:





验证的结论:
1. 13、21和18号染色体基因分布比较少,这三条染色体均可发生三体综合征


2. 基因分布与GC的分布好像还是有一些差别,老师ppt中的图不知道是怎么得到的...






本帖子中包含更多资源

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

x
回复 支持 5 反对 0

使用道具 举报

0

主题

20

帖子

124

积分

注册会员

Rank: 2

积分
124
发表于 6 天前 | 显示全部楼层
本帖最后由 y461650833y 于 2017-1-18 12:04 编辑

因为是做与猪相关研究,所以,针对Sus_scrofa.Sscrofa10.2.81.gtf进行处理
因为对单命令行perl不熟,所以只是照着运行了下。Homo_sapiens.GRCh38.87.chr.gtf文件的结果就不贴了。看下猪的:
cat Sus_scrofa.Sscrofa10.2.87.chr.gtf |perl -alne '{print if $F[2] eq "gene"}' |cut -f 1 |sort |uniq -c
   2209 1
    524 10
    422 11
   1145 12
   1534 13
   1370 14
    960 15
    444 16
    688 17
    512 18
   2147 2
   1444 3
   1254 4
   1172 5
   1941 6
   1651 7
    847 8
   1427 9
     37 MT
   1120 X
     13 Y

cat Sus_scrofa.Sscrofa10.2.87.chr.gtf |perl -alne '{next unless $F[2] eq "gene"; /gene_biotype "(.*?)";/; print $1}' |cut -f 1 |sort |uniq -c
     15 antisense
      2 IG_C_gene
      6 IG_J_gene
     15 IG_V_gene
      2 IG_V_pseudogene
     35 lincRNA
    768 miRNA
    171 misc_RNA
      2 Mt_rRNA
     22 Mt_tRNA
      5 non_coding
     80 processed_transcript
  19442 protein_coding
    555 pseudogene
    161 rRNA
    569 snoRNA
   1011 snRNA
自己编程:照着直播里面的稍微修改了一下 发现竟然可以运行!看来我抄的还是可以,至少没有抄错了!哈哈!!
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/perl -w
while(<>){
 chomp;
 @arr=split /\t/;
 if ($arr[2] eq "gene") {
 @brr = split  /\t/, $_;
 @crr = $brr[0];
 foreach $crr (@crr){
 $count{$crr}+=1;
 }
 }
}
foreach $crr(sort keys %count){
 print "$crr\t$count{$crr}\n";
}

各个染色体gene数量结果:
1       2209
10      524
11      422
12      1145
13      1534
14      1370
15      960
16      444
17      688
18      512
2       2147
3       1444
4       1254
5       1172
6       1941
7       1651
8       847
9       1427
MT      37
X       1120
Y       13
暂时只把这个做通了,剩下部分等学习明白了,在写!
吃完午饭,回来研究了一会代码,发现昨天写的代码,把{}放错地方了,真是尴尬!下面是看了直播的视频后,自己参考着写的!
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/perl -w
while(<>){
 chomp;
 @arr=split /\t/;
 if ($arr[2] eq "gene") {
  /gene_biotype "(.*?)";/;
 push @brr, $1;
 
 }
 }
 foreach $brr (@brr){
 $count{$brr}+=1;
 }
 

foreach $brr(sort keys %count){
 print "$brr\t$count{$brr}\n";
}

##gene_biotype结果:
IG_C_gene        2
IG_J_gene        6
IG_V_gene        15
IG_V_pseudogene        2
Mt_rRNA        2
Mt_tRNA        22
antisense        15
lincRNA        35
miRNA        768
misc_RNA        171
non_coding        5
processed_transcript        80
protein_coding        19442
pseudogene        555
rRNA        161
snRNA        1011
snoRNA        569
做完这个作业:发现自己要学习的知识还有很多,还要去慢慢学习了!特别是对于哈希的灵活运用,这个真是我的弱点,需要努力!
027-R+perl-游戏 一个perl小白,求解救!
人生若只如初见!
回复 支持 1 反对 0

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-12 22:37:20 | 显示全部楼层
[AppleScript] 纯文本查看 复制代码
import time
start = time.clock()
lichr = []
summa = {}
sum_gene = 0
sum_exon = 0
sum_transcript = 0

with open("F:\\Python34\\zbioinfo\\Homo_sapiens.GRCh38.87.chr.gtf","r") as f:
    for line in f:
        if line.startswith('#'):
            continue
        line = line.rstrip()
        chr_id = line.split('\t')[0]
        if chr_id not in lichr:
            lichr.append(chr_id)
            summa[chr_id] = {}
            summa[chr_id]["totitem"] = 0
            summa[chr_id]["gene"] = 0
            summa[chr_id]["exon"] = 0
            summa[chr_id]["transcript"] = 0
        summa[chr_id]["totitem"] += 1
        if "gene" == line.split('\t')[2]:
            summa[chr_id]["gene"] += 1
        if "exon" == line.split('\t')[2]:
            summa[chr_id]["exon"] += 1
        if "transcript" == line.split('\t')[2]:
            summa[chr_id]["transcript"] +=1

geneTotal = 0
exonTotal = 0
transcriptTotal = 0
for chr_id,sum_gene_exon_tran in summa.items():
    geneTotal += summa[chr_id]["gene"]
    exonTotal += summa[chr_id]["exon"]
    transcriptTotal += summa[chr_id]["transcript"]
    print (chr_id)
    print ("totitem : %s" % summa[chr_id]["totitem"])
    print ("gene : %s" % summa[chr_id]["gene"])
    print ("exon : %s" % summa[chr_id]["exon"])
    print ("transcript : %s" % summa[chr_id]["transcript"])
print ()    
print ("geneTotal :%s" % geneTotal)  
print ("exonTotal :%s" % exonTotal)
print ("transcriptTotal :%s" % transcriptTotal) 
end = time.clock()
print("The function run time is : %.03f seconds" %(end-start))
回复 支持 反对

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-12 22:39:04 | 显示全部楼层
上面的代码是python统计每个染色体中gene,exon,transcirpt的个数,感谢东哥的视频。
回复 支持 反对

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-12 22:40:54 | 显示全部楼层
[Python] 纯文本查看 复制代码
import time
start = time.clock()
lichr = []
ens_types = []
summ = {}
#i = 0
with open("F:\\Python34\\zbioinfo\\Homo_sapiens.GRCh38.87.chr.gtf","r") as f:
    for line in f:
        if line.startswith('#'):
            continue
        line = line.rstrip()
        chr_id = line.split('\t')[0]
        if chr_id not in lichr:
            lichr.append(chr_id)
        ens_type = line.split('\t')[2]
        if ens_type not in ens_types:
            ens_types.append(ens_type)
        
print (lichr)
print (ens_types)
end = time.clock()
print("The function run time is : %.03f seconds" %(end-start))


统计每个染色体中ensembl类型
回复 支持 反对

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-12 22:43:47 | 显示全部楼层
[AppleScript] 纯文本查看 复制代码
import sys
import time
start = time.clock()
lichr = []
ens_types = ['gene', 'transcript', 'exon', 'CDS', 'start_codon', 'stop_codon', 'five_prime_utr', 'three_prime_utr', 'Selenocysteine']
summ = {}
with open("F:\\Python34\\zbioinfo\\Homo_sapiens.GRCh38.87.chr.gtf","r") as f:
    for line in f:
        if line.startswith('#'):
            continue
        line = line.rstrip()
        chr_id = line.split('\t')[0]
        ens_type_chr = line.split('\t')[2]
        if chr_id not in lichr:
            lichr.append(chr_id) 
            summ[chr_id] = {}
            for ens_type in ens_types:
                summ[chr_id][ens_type] = 0
        for ens_type in ens_types:
            if ens_type_chr == ens_type:
                summ[chr_id][ens_type] += 1

for chr_id, sum_enstype in summ.items():
    print (chr_id)
    for ens_type, num in sum_enstype.items():
        print ("%s : %s" %(ens_type, num))
    print ()    
print ()
summary = {}
for ens_type in ens_types:
        summary[ens_type] = 0 
for chr_id,sum_enstype in summ.items():
    for ens_type,num in sum_enstype.items():
        summary[ens_type] += num
for ens_type, num in summary.items():
    print ('%s : %s' % (ens_type, num))
    
end = time.clock()
print("The function run time is : %.03f seconds" %(end-start))


统计每个染色体中ensembl类型的个数,统计所有染色体的ensembl类型的个数。谁能用python把统计染色体ensembl类型和类型的个数两个代码合到一起吗?我的能力不够。
回复 支持 反对

使用道具 举报

10

主题

32

帖子

139

积分

注册会员

Rank: 2

积分
139
发表于 2017-1-15 20:12:48 | 显示全部楼层
第三讲 Python-dongye代码

结果文件:

链接: http://pan.baidu.com/s/1qX89Pzm 密码: efev

代码:

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

args = sys.argv


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 = {}
    list_transcript = {}
    list_exon = []
    # l_n = 0
    with open(args[1]) as fp_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)
                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)


回复 支持 反对

使用道具 举报

8

主题

27

帖子

174

积分

注册会员

Rank: 2

积分
174
发表于 2017-1-15 20:26:07 | 显示全部楼层
本帖最后由 anlan 于 2017-1-22 18:15 编辑

034-perl-shell
基础部分-perl代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my %hash;
my %tmp;
open my $fh, "Homo_sapiens.GRCh38.87.chr.gtf" or die "Cannot open it";
while(<$fh>){
        chomp;
        my @array = split /\t/, $_;
        my $tmp;
        $hash{$array[0]} -> {$array[2]}++;
        if ($array[2] eq "gene"){
                if($array[8] =~ /gene_biotype.*?(protein_coding)/){
                        $hash{$array[0]} -> {$1}++; 
                }
                if($array[8] =~ /gene_biotype\s+\"(.*?)\"/){
                        $tmp{$1}++;
                }
        }
}
close $fh;

print "chr\tgene\tprotein_coding\ttranscript\texon\n";
map{
        print "$_\t$hash{$_}->{gene}\t$hash{$_}->{protein_coding}\t$hash{$_}->{transcript}\t$hash{$_}->{exon}\n";
}sort keys %hash;

print "\n\n";

map{
        print "$_\t$tmp{$_}\n";
}sort keys %tmp;

结果如下:
chr        gene        protein_coding        transcript        exon                        
1        5194        2052        17296        108417
10        2204        732        6507        41978
11        3235        1276        12151        69915
12        2940        1034        11419        69686
13        1304        327        3236        18266
14        2224        623        7476        41897
15        2152        597        7444        45517
16        2511        866        9896        57261
17        2995        1197        12573        74785
18        1170        270        3640        20309
19        2926        1470        12695        71424
2        3971        1255        14390        91110
20        1386        544        4242        25502
21        835        233        2413        13966
22        1318        438        4472        26063
3        3010        1076        12064        74562
4        2505        751        7732        45520
5        2868        876        9200        52715
6        2863        1045        8540        52060
7        2867        906        9437        56882
8        2353        676        7717        42730
9        2242        781        6581        42213
MT        37        13        37        37
X        2359        841        6037        35309
Y        523        53        740        3784


3prime_overlapping_ncRNA        30
IG_C_gene        14
IG_C_pseudogene        9
IG_D_gene        37
IG_J_gene        18
IG_J_pseudogene        3
IG_V_gene        145
IG_V_pseudogene        193
IG_pseudogene        1
Mt_rRNA        2
Mt_tRNA        22
TEC        1048
TR_C_gene        6
TR_D_gene        4
TR_J_gene        79
TR_J_pseudogene        4
TR_V_gene        108
TR_V_pseudogene        30
antisense        5526
bidirectional_promoter_lncRNA        4
lincRNA        7533
macro_lncRNA        1
miRNA        1567
misc_RNA        2212
non_coding        3
polymorphic_pseudogene        51
processed_pseudogene        10268
processed_transcript        515
protein_coding        19932
pseudogene        21
rRNA        543
ribozyme        8
sRNA        5
scRNA        1
scaRNA        49
sense_intronic        908
sense_overlapping        187
snRNA        1900
snoRNA        944
transcribed_processed_pseudogene        450
transcribed_unitary_pseudogene        60
transcribed_unprocessed_pseudogene        735
unitary_pseudogene        154
unprocessed_pseudogene        2661
vaultRNA        1

高级部分代码如下:
gene_trans函数为基因to转录本的分布统计函数
tran_exons函数为转录本to外显子的分布统计函数

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my @in = ("Human.gtf", "Mus.gtf", "Rat.gtf", "Dog.gtf", "Cat.gtf");
foreach my $file (@in){
        my @file_split = split /\./, $file;
        my $file_out = $file_split[0].".txt";
        open my $fh_out, ">$file_out" or die "Cannot write it";
        my %tmp = &tran_exons($file);
        my $all;
        map{
                $all += $tmp{$_} if $_ >0 
        }keys %tmp;
        map{
                my $percent = $tmp{$_}/$all;
                print $fh_out "$_\t$tmp{$_}\t$percent\n" if $_ > 0;
        }sort{$a<=>$b} keys %tmp;
        close $fh_out;
}

sub gene_trans{
        my $file_trans = shift @_;
        open my $fh, $file_trans or die "Cannot open $file_trans";
        my $i = 0;
        my %hash;
        while(<$fh>){
                my @array = split /\t/, $_;
                next if $_ =~ /^\#/;
                if ($array[2] && $array[2] eq "gene"){
                        if ($hash{$i}){
                                $hash{$i}++;
                        }else{
                                $hash{$i} = 1;
                        }
                        $i = 0;
                }else{
                        if ($array[2] eq "transcript"){
                                $i++;
                        }
                }
        }
        close $fh;
        return %hash;
}


sub tran_exons{
        my $file_trans = shift @_;
        open my $fh, $file_trans or die "Cannot open $file_trans";
        my $i = 0;
        my %hash;
        while(<$fh>){
                my @array = split /\t/, $_;
                next if $_ =~ /^\#/;
                if ($array[2] && $array[2] eq "transcript"){
                        if ($hash{$i}){
                                $hash{$i}++;
                        }else{
                                $hash{$i} = 1;
                        }
                        $i = 0;
                }else{
                        if ($array[2] eq "exon"){
                                $i++;
                        }
                }
        }
        close $fh;
        return %hash;
}


shell部分代码:每个染色体上gene数目统计:
[Shell] 纯文本查看 复制代码
cat Human.gtf | grep -v '^#' |awk '$3~/gene/{print $1,"\t",$3}' | sort |uniq -c |less -S

每个染色体上protein coding gene统计:
[Shell] 纯文本查看 复制代码
cat Human.gtf | grep -v '^#' |awk '$3~/gene/{print $0}' | awk '$0~/protein_coding/{print $1,"\t", $3}' | sort | uniq -c |less -S
每个基因的类型统计:
还是不确定awk 的 match函数有无非贪婪匹配,等以后正式学shell再来修改
错误的代码:
[Shell] 纯文本查看 复制代码
cat Human.gtf | grep -v '^#' |awk '$3~/gene/{print $0}' | awk 'match($0, /gene_biotype\s+"(.+?)";/, a) {print a[1]}'|less -S 

正确的代码:
[Shell] 纯文本查看 复制代码
cat Human.gtf | grep -v '^#' |awk '$3~/gene/{print $0}' | awk 'match($0, /gene_biotype\s+"(.+?)";(\shavana_gene\s)/, a) {print a[1]}'| sort | uniq -c |less -S

第二个代码是把匹配更加精确了,减少了第一个代码引起的贪婪匹配

回复 支持 反对

使用道具 举报

0

主题

8

帖子

59

积分

注册会员

Rank: 2

积分
59
发表于 7 天前 | 显示全部楼层
本帖最后由 xxnku 于 2017-1-17 11:38 编辑

简单的统计gene numbers
[Python] 纯文本查看 复制代码
import re
import time
from collections import OrderedDict
start = time.clock()


chr_id = ''
tmp_dict = OrderedDict()
with open('/home/zijun/biotrainee/third/Homo_sapiens.GRCh38.87.chr.gtf','r') as f:
    for line in f.readlines():
        if line.startswith('#'):
            continue
        line = line.rstrip().split('\t')
        chr_id = line[0]

        if chr_id not in tmp_dict.keys():
            tmp_dict[chr_id] = 0

        if re.match('gene',line[2]):
            tmp_dict[chr_id] += 1

end = time.clock()
print ("%.5f"%(end-start))

for k,v in tmp_dict.items():
    print (k,v)


感觉自己写的代码好low...
回复 支持 反对

使用道具 举报

0

主题

4

帖子

38

积分

新手上路

Rank: 1

积分
38
发表于 7 天前 | 显示全部楼层
137+python
题目:简单输出染色体中gene、transcript、exon、cds的数量。
状态:看了视频,根据视频修改了思路,思考后结果。
思路:
数据结构:{chr_id:{feature_id:number}}
1 下载数据:wget
2 输入方法:sys.argv[1]
3 提取每一列至list[]
4 list[0]→chr_id;list[2]→feature_id
5 根据feature_id相加,得到结果
6 输出:open,也可以用argv[2]


[Python] 纯文本查看 复制代码
import time
import sys
from collections import OrderedDict
time_start = time.clock()

chr_gene = OrderedDict()
list_chr = []
list_feature = ['gene', 'transcript', 'exon', 'CDS']

with open(sys.argv[1], 'r') as file:
    for line in file:
        if line.startswith('#'):
            continue
        line = line.strip("\n")
        line = line.split("\t")
        chr_id = line[0]
        feature_id = line[2]
        if not chr_id in list_chr:
            list_chr.append(chr_id)
            chr_gene[chr_id] = {}
            for fea in list_feature:
                chr_gene[chr_id][fea] = 0
        for feature_id in list_feature:
            chr_gene[chr_id][feature_id] += 1

sorted(chr_gene.keys())

with open('chr_geneNum.txt',"w") as output:
    output.write("chr_item" + "\t" + "gene_num" + "\n")
    for chr_item, gene_num in chr_gene.items():
        chr_item = str(chr_item)
        gene_num = str(gene_num)
        # print("\t".join([chr_item, gene_num]) + "\n")
        output.write("\t".join([chr_item, gene_num]) + "\n")

time_end = time.clock()
print("Time : %.2f seconds" % (time_end - time_start))

结果:
chr_item        gene_num
1        {'CDS': 231832, 'exon': 231832, 'gene': 231832, 'transcript': 231832}
2        {'CDS': 192181, 'exon': 192181, 'gene': 192181, 'transcript': 192181}
3        {'CDS': 159878, 'exon': 159878, 'gene': 159878, 'transcript': 159878}
4        {'CDS': 99634, 'exon': 99634, 'gene': 99634, 'transcript': 99634}
5        {'CDS': 114524, 'exon': 114524, 'gene': 114524, 'transcript': 114524}
6        {'CDS': 115866, 'exon': 115866, 'gene': 115866, 'transcript': 115866}
7        {'CDS': 122010, 'exon': 122010, 'gene': 122010, 'transcript': 122010}
X        {'CDS': 80423, 'exon': 80423, 'gene': 80423, 'transcript': 80423}
8        {'CDS': 94090, 'exon': 94090, 'gene': 94090, 'transcript': 94090}
9        {'CDS': 90778, 'exon': 90778, 'gene': 90778, 'transcript': 90778}
11        {'CDS': 155702, 'exon': 155702, 'gene': 155702, 'transcript': 155702}
10        {'CDS': 90507, 'exon': 90507, 'gene': 90507, 'transcript': 90507}
12        {'CDS': 153049, 'exon': 153049, 'gene': 153049, 'transcript': 153049}
13        {'CDS': 38501, 'exon': 38501, 'gene': 38501, 'transcript': 38501}
14        {'CDS': 92999, 'exon': 92999, 'gene': 92999, 'transcript': 92999}
15        {'CDS': 97141, 'exon': 97141, 'gene': 97141, 'transcript': 97141}
16        {'CDS': 124801, 'exon': 124801, 'gene': 124801, 'transcript': 124801}
17        {'CDS': 164830, 'exon': 164830, 'gene': 164830, 'transcript': 164830}
18        {'CDS': 44609, 'exon': 44609, 'gene': 44609, 'transcript': 44609}
20        {'CDS': 55735, 'exon': 55735, 'gene': 55735, 'transcript': 55735}
19        {'CDS': 163169, 'exon': 163169, 'gene': 163169, 'transcript': 163169}
Y        {'CDS': 7282, 'exon': 7282, 'gene': 7282, 'transcript': 7282}
22        {'CDS': 56203, 'exon': 56203, 'gene': 56203, 'transcript': 56203}
21        {'CDS': 28892, 'exon': 28892, 'gene': 28892, 'transcript': 28892}
MT        {'CDS': 144, 'exon': 144, 'gene': 144, 'transcript': 144}


小结:
1 巩固open 的输入输出方法
2 尝试将目标转化成思路,并完成编程
3 简单数据结构整理
4 字典应用
问题:
不会将字典中嵌套的内容也按照表格的形式输出,待学习。

下一步:
完成python教学视频内容学习并写一遍,初步查看,内容包括:
1 def / class
2 git软件
3 复杂数据结构整理




回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-1-24 10:56 , Processed in 0.267192 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.