搜索
楼主: Jimmy

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

  [复制链接]

0

主题

5

帖子

115

积分

注册会员

Rank: 2

积分
115
发表于 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";
}

                  
                  
回复 支持 反对

使用道具 举报

0

主题

9

帖子

79

积分

注册会员

Rank: 2

积分
79
发表于 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的个数,看着老师的视频还有很多不太明白,仍需努力学习
回复 支持 反对

使用道具 举报

0

主题

9

帖子

79

积分

注册会员

Rank: 2

积分
79
发表于 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的个数,看着老师的视频还有很多不太明白,仍需努力学习
回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 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()
回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 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 ...

写成第四题了我去。。。
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 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等命令的基本用法
回复 支持 反对

使用道具 举报

13

主题

33

帖子

243

积分

版主

Rank: 7Rank: 7Rank: 7

积分
243
发表于 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))))


回复 支持 反对

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 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))
回复 支持 反对

使用道具 举报

0

主题

6

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 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作图
回复 支持 反对

使用道具 举报

1

主题

17

帖子

268

积分

中级会员

Rank: 3Rank: 3

积分
268
发表于 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


回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-21 01:17 , Processed in 0.044757 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.