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

033-Perl-R-Pyth · 发表于 2017-1-26 17:46:20

本帖最后由 033-Perl-R-Pyth 于 2017-1-26 17:47 编辑

看视频前
思路：
用hash区别不同的chr，用split函数选取目标列，遍历hash输出

[Perl] 纯文本查看 复制代码

#!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";
}

haozhou · 发表于 2017-1-27 16:26:15

010-python

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

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)

统计了一下每个染色体上gene，exon，transcript的个数，看着老师的视频还有很多不太明白，仍需努力学习

haozhou · 发表于 2017-1-27 16:26:42

010-python

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

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)

统计了一下每个染色体上gene，exon，transcript的个数，看着老师的视频还有很多不太明白，仍需努力学习

017中大普外 · 发表于 2017-1-28 15:14:55

写了一个TCGA版本的，结合了TCGA的metadata

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

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()

017中大普外 · 发表于 2017-1-28 20:57:14

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

写成第四题了我去。。。

熊海清 · 发表于 2017-1-30 23:28:23

本帖最后由熊海清于 2017-1-30 23:36 编辑

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

Andyhigh · 发表于 2017-2-1 16:01:16

092-R-python

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

[AppleScript] 纯文本查看 复制代码

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))))

Teanna2017 · 发表于 2017-2-4 16:56:29

本帖最后由 Teanna2017 于 2017-2-4 16:58 编辑

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

代码如下：

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

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))

雨林木风 · 发表于 2017-2-5 18:24:28

python-031
思路：
1. 按行读取数据，跳过#开头行
2. 提取每一行的信息至列表l_list
3. 构建数据结构：{chr_id：{feature:number}};
4. 输出统计结果
代码：

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

#! 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][i] += 1
                if i not in gtfcount[Chr]:
                    gtfcount[Chr][i] = 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"

结果：
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 · 发表于 2017-2-5 21:45:02

Python-048
先简单的统计gtf中每条染色体的gene数量，transcript数量，exon数量
思路:构建一个字典，以染色体号位key，gene数量，transcript数量，exon数量组成的列表为value

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

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

在输出列表时，想以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

		自动登录	找回密码
密码			立即注册