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

033-Perl-R-Pyth

看视频前
思路：
用hash区别不同的chr，用split函数选取目标列，遍历hash输出[mw_shl_code=perl,true]
#!usr/bin/perl -w
# GTF caculation
use strict; use warnings;
open A, 'Homo_sapiens.GRCh38.87.chr.gtf'or die $!;
my ($chr, %gene_count,%exon_count);
while (<A>){
        chomp;
        if (/^#/){
                next;
        }
        else{
                        $chr=(split(/\t/,$_))[0];
                    if( (split("\t",$_))[2] eq 'gene'){
                       $gene_count{$chr}=$gene_count{$chr}+1;
        }
                elsif ( (split("\t",$_))[2] eq 'exon'){
                        $exon_count{$chr}=$exon_count{$chr}+1;
                    }
}
                    }                  
close A;

foreach (sort keys %gene_count){
        #print "$_\n";
        print "chr$_\tgene\t$gene_count{$_}\n";
        print "chr$_\texon\t$exon_count{$_}\n";
}[/mw_shl_code]
                  
                  

haozhou

010-python
[mw_shl_code=python,true]import re
import os
from collections import OrderedDict
from operator import itemgetter
os.chdir('C:/Users/Administrator/Desktop/')

genesum=dict()
exonsum=dict()
transcriptsum=dict()
with open ('100hg19gtf') as f:
    for line in f:
        if line.startswith('#'):
            continue
        line=line.strip()
        #print(line)
        word=line.split()
        if word[2]=='gene':
            #if word[0] not in genesum:
                #genesum[word[0]]=1
            #else:
                #genesum[word[0]]+=1
#这里的if。。else 和下面的get作用一样的，get简洁一些
            genesum[word[0]]=genesum.get(word[0],0)+1
        elif word[2]=='exon':
            exonsum[word[0]]=exonsum.get(word[0],0)+1
        elif word[2]=='transcript':
            transcriptsum[word[0]]=transcriptsum.get(word[0],0)+1
        
print (genesum)
print(exonsum)
print(transcriptsum)       [/mw_shl_code]统计了一下每个染色体上gene，exon，transcript的个数，看着老师的视频还有很多不太明白，仍需努力学习

haozhou

010-python
[mw_shl_code=python,true]import re
import os
from collections import OrderedDict
from operator import itemgetter
os.chdir('C:/Users/Administrator/Desktop/')

genesum=dict()
exonsum=dict()
transcriptsum=dict()
with open ('100hg19gtf') as f:
    for line in f:
        if line.startswith('#'):
            continue
        line=line.strip()
        #print(line)
        word=line.split()
        if word[2]=='gene':
            #if word[0] not in genesum:
                #genesum[word[0]]=1
            #else:
                #genesum[word[0]]+=1
#这里的if。。else 和下面的get作用一样的，get简洁一些
            genesum[word[0]]=genesum.get(word[0],0)+1
        elif word[2]=='exon':
            exonsum[word[0]]=exonsum.get(word[0],0)+1
        elif word[2]=='transcript':
            transcriptsum[word[0]]=transcriptsum.get(word[0],0)+1
        
print (genesum)
print(exonsum)
print(transcriptsum)       [/mw_shl_code]统计了一下每个染色体上gene，exon，transcript的个数，看着老师的视频还有很多不太明白，仍需努力学习

017中大普外

写了一个TCGA版本的，结合了TCGA的metadata[mw_shl_code=python,true]import os
import sys
from collections import OrderedDict
import json
import re
import glob


gene_group=OrderedDict()
tcga_file=OrderedDict()
list_1=[]
normal=[]
tumor=[]
os.chdir("E:/z/files")

#处理TCGA metadata
with open('metadata.cart.2016-11-03T07_53_29.306026.json',"r") as f1:
    data=json.load(f1)
    for i in data:
        file_name=i["file_name"]
        file_name=re.sub("\.gz","",file_name)
        tcga_id=i["associated_entities"][0]["entity_submitter_id"]
        tcga_file[tcga_id]=file_name
        ID=tcga_id.split("-")
        if ID[3]=="11A":
            normal.append(tcga_id)              
        else:
            tumor.append(tcga_id)
list_1.extend(normal)
list_1.extend(tumor)
list_2=[]
for items in list_1:
    tcgafile=tcga_file[items]
    list_2.append(tcgafile)
##y=list_2[1]
#print(x,y)
def readfile(file):
    file_obj=open(file)
    try:
        while True:
            line=file_obj.readline().strip("\n")
            if not line:
                break
            arry=line.split("\t")
            if arry[0] in gene_group:
                gene_group[arry[0]].append(arry[1])
            else:
                gene_group[arry[0]]=[arry[1]]
    finally:
        file_obj.close()
    return gene_group
f3=open("TCGA-rnaseq-matrix.txt","w")
def main():
    for file in list_2:
        if glob.glob(file):           
            readfile(file)
    for gene_name,gene_expression in gene_group.items():
        
        f3.write(gene_name+"\t"+"\t".join(gene_expression)+'\n')
        
if __name__=="__main__":
    main()[/mw_shl_code]

017中大普外

017中大普外 发表于 2017-1-28 15:14
写了一个TCGA版本的，结合了TCGA的metadata[mw_shl_code=python,true]import os
import sys
from collectio ...

写成第四题了我去。。。

熊海清

074-python
1.本次视频中老师的讲解和他用的代码没有完全理解，再学习一段时间再补做本次python练习
2.学习了解了人类基因组GTF注释文件的内容
3.用shell练习了这次的习题，统计了人类基因组基因的数目和蛋白编码基因的数目。熟悉了shell中cat,grep,awk正则表达式,sort,uniq等命令的基本用法

Andyhigh

092-R-python

自己学习了下 用R解决问题

[mw_shl_code=applescript,true]rm(list=ls())
setwd("C:/Users/Administrator/Desktop/learn/")

source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)

###先赋值，再统计  转录本
x <- org.Hs.egENSEMBLTRANS

# Get the entrez gene IDs that are mapped to an Ensembl ID
mapped_genes <- mappedkeys(x)

# Convert to a list
xx <- as.list(x[mapped_genes])

length(x)
###可知共有 个基因都是有entrez ID号能对应转录本ID的基因有多少个呢？

length(xx)
#可以看到共有8055个基因都是有转录本的！？？？


###能对应ensembl的gene ID的基因有多少个呢？
a <- org.Hs.egENSEMBL

# Get the entrez gene IDs that are mapped to an Ensembl ID
mapped_genes <- mappedkeys(a)

# Convert to a list
aa <- as.list(a[mapped_genes])

length(a)
length(aa)

##有25725个基因 有对应的ensembl中gene ID  
###那么基因对应的转录本分布情况如何呢？

table(unlist(lapply(aa,length)))

###可以看出绝大部分的基因都是20个转录本一下的，但也有极个别基因居然有高达两百个转录本，很可怕！
##那么基因在染色体的分布情况如何呢？



b <- org.Hs.egCHR

# Get the entrez gene identifiers that are mapped to a chromosome

mapped_genes <- mappedkeys(b)

# Convert to a list

bb <- as.list(b[mapped_genes])

length(b)
length(bb)

###可以看到有接近五百个基因居然是没有染色体定位信息的！！！


table(unlist(bb))
bbb<-table(unlist(bb))
barplot(bbb)
##用barplot函数可视化一下，如图


###那么有多多少基因是有GO注释的呢？
c <- org.Hs.egGO

# Get the entrez gene identifiers that are mapped to a GO ID
mapped_genes <- mappedkeys(c)

# Convert to a list
cc <- as.list(c[mapped_genes])

length(cc)
length(c)

##可以看到只有18622个基因是有go注释信息的。
##那么基因被注释的go的分布如何呢？

table(unlist(lapply(cc,length)))
###可以看到大部分的基因都是只有30个go的，
##但是某些基因特别活跃，高达197个go注释。
barplot(table(unlist(lapply(cc,length))))
[/mw_shl_code]

Teanna2017

python-60

对每条染色体上的基因，转录本及外显子个数，平均每个基因有几个转录本，平均每个转录本有几个外显子进行了统计，脚本中还对所有染色体的基因，转录本，外显子，及平均每个基因的转录本数，平均转录本的外显子数进行了统计，发现1号染色体的基因数目最多，而Y染色体的基因数目最少，具体结果如下：

chr_1 gene 5194
chr_1 transcript/gene 3.33
chr_1 exon/transcript 6.27
chr_2 gene 3971
chr_2 transcript/gene 3.62
chr_2 exon/transcript 6.33
chr_3 gene 3010
chr_3 transcript/gene 4.01
chr_3 exon/transcript 6.18
chr_4 gene 2505
chr_4 transcript/gene 3.09
chr_4 exon/transcript 5.89
chr_5 gene 2868
chr_5 transcript/gene 3.21
chr_5 exon/transcript 5.73
chr_6 gene 2863
chr_6 transcript/gene 2.98
chr_6 exon/transcript 6.1
chr_7 gene 2867
chr_7 transcript/gene 3.29
chr_7 exon/transcript 6.03
chr_X gene 2359
chr_X transcript/gene 2.56
chr_X exon/transcript 5.85
chr_8 gene 2353
chr_8 transcript/gene 3.28
chr_8 exon/transcript 5.54
chr_9 gene 2242
chr_9 transcript/gene 2.94
chr_9 exon/transcript 6.41
chr_11 gene 3235
chr_11 transcript/gene 3.76
chr_11 exon/transcript 5.75
chr_10 gene 2204
chr_10 transcript/gene 2.95
chr_10 exon/transcript 6.45
chr_12 gene 2940
chr_12 transcript/gene 3.88
chr_12 exon/transcript 6.1
chr_13 gene 1304
chr_13 transcript/gene 2.48
chr_13 exon/transcript 5.64
chr_14 gene 2224
chr_14 transcript/gene 3.36
chr_14 exon/transcript 5.6
chr_15 gene 2152
chr_15 transcript/gene 3.46
chr_15 exon/transcript 6.11
chr_16 gene 2511
chr_16 transcript/gene 3.94
chr_16 exon/transcript 5.79
chr_17 gene 2995
chr_17 transcript/gene 4.2
chr_17 exon/transcript 5.95
chr_18 gene 1170
chr_18 transcript/gene 3.11
chr_18 exon/transcript 5.58
chr_20 gene 1386
chr_20 transcript/gene 3.06
chr_20 exon/transcript 6.01
chr_19 gene 2926
chr_19 transcript/gene 4.34
chr_19 exon/transcript 5.63
chr_Y gene 523
chr_Y transcript/gene 1.41
chr_Y exon/transcript 5.11
chr_22 gene 1318
chr_22 transcript/gene 3.39
chr_22 exon/transcript 5.83
chr_21 gene 835
chr_21 transcript/gene 2.89
chr_21 exon/transcript 5.79
chr_MT gene 37
chr_MT transcript/gene 1.0
chr_MT exon/transcript 1.0
all_chr_gene 57992
all_chr_transcript 197935
all_exon 1181908
all_chr_transcript/gene 3.41
all_chr_exon/transcript 5.97

代码如下：

[mw_shl_code=python,true]import re
import collections

path = '/Users/sunshine/practice/practice3-hg38/'

f = open(path+'Homo_sapiens.GRCh38.87.chr.gtf','r')
f_lines = f.readlines()

chr_dict = collections.OrderedDict()  #建立有序字典
temps = []

for line in f_lines:
    line = line.strip('\n')
    if line.startswith('#'):
        pass
    else:
        line = line.split('\t')
        chr_id = 'chr_'+line[0]
        chr_dict[chr_id] = {}
        temps.append(line[2])  #对文件中第三列的每种类型进行合并
        temps=list(set(temps))  #合并，去重
        for temp in temps:
            chr_dict[chr_id][temp] = 0
        
for line in f_lines:
    line = line.strip('\n')
    if line.startswith('#'):
        pass
    else:
        line = line.split('\t')
        chr_id = 'chr_'+line[0]
        chr_dict[chr_id][line[2]] += 1 #这里的line[2]即temp，故这里为对每个染色体中的每种情况进行个数统计
        
gene_sum = 0
transcript_sum = 0
exon_sum = 0

for chr_id, situations in chr_dict.items():
    gene_sum  += situations['gene']
    transcript_sum += situations['transcript']
    exon_sum += situations['exon']
    print (chr_id,'gene',situations['gene'])
    print (chr_id,'transcript/gene',round(situations['transcript']/situations['gene'],2))
    print (chr_id,'exon/transcript',round(situations['exon']/situations['transcript'],2))

print ('all_chr_gene',gene_sum)
print ('all_chr_transcript',transcript_sum)
print ('all_exon',exon_sum)
print ('all_chr_transcript/gene',round(transcript_sum/gene_sum,2))
print ('all_chr_exon/transcript',round(exon_sum/transcript_sum,2))
[/mw_shl_code]

雨林木风

python-031
思路：
1. 按行读取数据，跳过#开头行
2. 提取每一行的信息至列表l_list
3. 构建数据结构：{chr_id：{feature:number}};
4. 输出统计结果
代码：
[mw_shl_code=python,true]#! usr/bin/envpython
# _*_coding:utf-8_*_

from collections import OrderedDict
import sys
import time

start = time.clock()
gtfcount = OrderedDict()
featurelist = ['gene','transcript','CDS','exon']
with open(sys.argv[1], 'r') as f:
    for line in f:
        if not line.startswith('#'):   
            line = line.strip()
            l_list = line.split("\t")
            Chr = l_list[0]
            feature = l_list[2]
            if Chr not in gtfcount.keys():
                gtfcount[Chr] = {}
            for i in featurelist:
                if not feature == i:
                    continue
                if i in gtfcount[Chr]:
                    gtfcount[Chr] += 1
                if i not in gtfcount[Chr]:
                    gtfcount[Chr] = 1
with open("result.txt", 'w') as fp_out:
    fp_out.write("\t".join(["Chr","gtfcount"]) + "\n")
    print "\t".join(["Chr","gtfcount"]) + "\n"
    for Chr, gtfcount in gtfcount.items():
        fp_out.write("\t".join([Chr, str(gtfcount)]) + "\n")
        print "\t".join([Chr,str(gtfcount)]) + "\n"
[/mw_shl_code]

结果：
Chr     gtfcount
1       {'exon': 108417, 'gene': 5194, 'transcript': 17296, 'CDS': 64975}
2       {'exon': 91110, 'gene': 3971, 'transcript': 14390, 'CDS': 53480}
3       {'exon': 74562, 'gene': 3010, 'transcript': 12064, 'CDS': 43036}
4       {'exon': 45520, 'gene': 2505, 'transcript': 7732, 'CDS': 27178}
5       {'exon': 52715, 'gene': 2868, 'transcript': 9200, 'CDS': 30758}
6       {'exon': 52060, 'gene': 2863, 'transcript': 8540, 'CDS': 33020}
7       {'exon': 56882, 'gene': 2867, 'transcript': 9437, 'CDS': 32904}
X       {'exon': 35309, 'gene': 2359, 'transcript': 6037, 'CDS': 22767}
8       {'exon': 42730, 'gene': 2353, 'transcript': 7717, 'CDS': 24440}
9       {'exon': 42213, 'gene': 2242, 'transcript': 6581, 'CDS': 26270}
11      {'exon': 69915, 'gene': 3235, 'transcript': 12151, 'CDS': 42363}
10      {'exon': 41978, 'gene': 2204, 'transcript': 6507, 'CDS': 26710}
12      {'exon': 69686, 'gene': 2940, 'transcript': 11419, 'CDS': 42166}
13      {'exon': 18266, 'gene': 1304, 'transcript': 3236, 'CDS': 10435}
14      {'exon': 41897, 'gene': 2224, 'transcript': 7476, 'CDS': 24528}
15      {'exon': 45517, 'gene': 2152, 'transcript': 7444, 'CDS': 26068}
16      {'exon': 57261, 'gene': 2511, 'transcript': 9896, 'CDS': 32673}
17      {'exon': 74785, 'gene': 2995, 'transcript': 12573, 'CDS': 44869}
18      {'exon': 20309, 'gene': 1170, 'transcript': 3640, 'CDS': 11850}
20      {'exon': 25502, 'gene': 1386, 'transcript': 4242, 'CDS': 15499}
19      {'exon': 71424, 'gene': 2926, 'transcript': 12695, 'CDS': 43907}
Y       {'exon': 3784, 'gene': 523, 'transcript': 740, 'CDS': 1538}
22      {'exon': 26063, 'gene': 1318, 'transcript': 4472, 'CDS': 14888}
21      {'exon': 13966, 'gene': 835, 'transcript': 2413, 'CDS': 7396}
MT      {'exon': 37, 'gene': 37, 'transcript': 37, 'CDS': 13}

小结：看了东哥的视频很受启发，了解了写一个可复用的程序的思路，也从这次视频中熟悉了面向对象编程，虽然自己目前还写不出来这样的程序，但是仍然给了我今后处理问题的很好的思路。
下一步：学习怎样用python及R作图

btrainee

Python-048
先简单的统计gtf中每条染色体的gene数量，transcript数量，exon数量
思路:构建一个字典，以染色体号位key，gene数量，transcript数量，exon数量组成的列表为value
[mw_shl_code=python,true]#! /usr/bin/env python3
#! -*- coding:utf-8 -*-

import sys
from collections import OrderedDict

args = sys.argv

def main(args):
        CHR_dict = OrderedDict()
        with open(args[1],'r') as fp_gtf:
                for line in fp_gtf:
                        if line.startswith('#'):
                                continue
                        lines = line.strip().split('\t')
                        chr = lines[0]
                        if chr not in CHR_dict:
                                CHR_dict[chr] = [0,0,0]
                        elif lines[2] == "gene":
                                CHR_dict[chr][0] += 1
                        elif lines[2] == "transcript":
                                CHR_dict[chr][1] += 1
                        elif lines[2] == "exon":
                                CHR_dict[chr][2] += 1

        with open(args[2],'w') as fp_out:
                fp_out.write("Chr\tGene_num\tTranscript_num\tExon_num\n")
                for chr,info in CHR_dict.items():
                        fp_out.write("\t".join([chr,str(info[0]),str(info[1]),str(info[2])]) + "\n")

if __name__ == "__main__":
        main(args)[/mw_shl_code]
在输出列表时，想以tab键分割，目前只想到这种解决方式"\t".join([chr,str(info[0]),str(info[1]),str(info[2])]) 。当列表值比较多时就比较费劲，
应该还会有更简单的方式
基因组（Homo_sapiens.GRCh38.86.gtf），部分结果如下：
Chr    Gene_num    Transcript_num    Exon_num
1    5193    17296    108417
2    3970    14390    91110
3    3009    12064    74562
4    2504    7732    45520
5    2867    9200    52715
6    2862    8540    52060
7    2866    9437    56882
X    2358    6037    35309
8    2352    7717    42730
9    2241    6581    42213
11    3234    12151    69915
10    2203    6507    41978
12    2939    11419    69686
13    1303    3236    18266
14    2223    7476    41897
15    2151    7444    45517