搜索
楼主: Jimmy

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

  [复制链接]

0

主题

2

帖子

40

积分

新手上路

Rank: 1

积分
40
发表于 2017-2-14 14:46:50 | 显示全部楼层

128-R+perl+shell-嗖

本帖最后由 ynsu 于 2017-2-14 14:48 编辑

感慨还是挺多的 ,不敢轻易写,看了shell,perl视频,又好好思考了下,彻底搞透彻编程思想才敢下手写的,编程思想真的特别特别重要,更进一步体会到了Jimmy及群里的前辈们说的话了,编程思想懂了,多种语言均可实现编程的目的,从小伙伴儿们的作业中也能学到自己所不知道的新知识,再次感谢,谢谢优秀的平台和优秀的小伙伴们!
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my %hashone;my %hashtwo;my @t;
open IN,$ARGV[0];
while (<IN>) {
	chomp;
	next if /^#/;
	@t=split/\t/,$_;
	$hashone{$t[0]}++ if $t[2]=~/gene/;	
	$hashtwo{$t[0]}++ if $t[2]=~/transcript/ && $_=~/protein_coding/;
	}
close IN;
	foreach my $chr(sort keys %hashone){
		print "$chr\tgene\t$hashone{$chr}\tprotein_coding\t$hashtwo{$chr}\n";
	}
 


附下结果:
1       gene    5194    protein_coding  12869
10      gene    2204    protein_coding  4421
11      gene    3235    protein_coding  9633
12      gene    2940    protein_coding  8925
13      gene    1304    protein_coding  1868
14      gene    2224    protein_coding  5446
15      gene    2152    protein_coding  5258
16      gene    2511    protein_coding  7707
17      gene    2995    protein_coding  10020
18      gene    1170    protein_coding  2453
19      gene    2926    protein_coding  10680
2       gene    3971    protein_coding  9955
20      gene    1386    protein_coding  3084
21      gene    835     protein_coding  1497
22      gene    1318    protein_coding  3213
3       gene    3010    protein_coding  8932
4       gene    2505    protein_coding  5424
5       gene    2868    protein_coding  6576
6       gene    2863    protein_coding  6071
7       gene    2867    protein_coding  6739
8       gene    2353    protein_coding  5497
9       gene    2242    protein_coding  4652
MT      gene    37      protein_coding  13
X       gene    2359    protein_coding  4258
Y       gene    523     protein_coding  195

实在不好意思!总发错地方!!!以后一定注意!忘大家谅解!
回复 支持 反对

使用道具 举报

2

主题

21

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-2-15 10:51:12 | 显示全部楼层
本帖最后由 AdaWon 于 2017-2-15 11:00 编辑

147-python
#本次任务是:gff3/gtf注释文件的探究
#选用 Arabidopsis 注释文件
1.查看注释文件;
$ less TAIR10_GFF3_genes.gff
Chr1    TAIR10  chromosome      1       30427671        .       .       .       ID=Chr1;Name=Chr1
Chr1    TAIR10  gene    3631    5899    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1    TAIR10  mRNA    3631    5899    .       +       .       ID=AT1G01010.1; Parent=AT1G01010;Name=AT1G01010.1;Index=1
Chr1    TAIR10  protein 3760    5630    .       +       .       ID=AT1G01010.1-Protein;Name=AT1G01010.1;Derives_from=AT1G01010.1
Chr1    TAIR10  exon    3631    3913    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  five_prime_UTR  3631    3759    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     3760    3913    .       +       0       Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1    TAIR10  exon    3996    4276    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     3996    4276    .       +       2       Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1    TAIR10  exon    4486    4605    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     4486    4605    .       +       0       Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1    TAIR10  exon    4706    5095    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     4706    5095    .       +       0       Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1    TAIR10  exon    5174    5326    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     5174    5326    .       +       0       Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1    TAIR10  exon    5439    5899    .       +       .       Parent=AT1G01010.1
Chr1    TAIR10  CDS     5439    5630    .       +       0       Parent=AT1G01010.1,AT1G01010.1-Protein;


$ cut -f3 TAIR10_GFF3_genes.gff |sort |uniq
CDS
chromosome
exon
five_prime_UTR
gene
miRNA
mRNA
mRNA_TE_gene
ncRNA
protein
pseudogene
pseudogenic_exon
pseudogenic_transcript
rRNA
snoRNA
snRNA
three_prime_UTR
transposable_element_gene
tRNA


2. 根据注释文件的内容,对python教学视频的代码略作调整,以实现对拟南芥注释文件的探索;

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


import sys
import re #在Python中使用正则表达式需要标准库中的一个包re

args = sys.argv 
#sys.argv[]是用来获取命令行参数的,sys.argv[0]表示代码本身文件路径,所以参数从1开始

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 = []

        with open(args[1]) as fp_gtf:
                for line in fp_gtf:
                        if line.startswith("#"): #注意:是starts 非 start
                                #startswith() 方法用于检查字符串是否是以指定子字符串开头,如果是则返回 True,否则返回 False。
                                continue
                                #continue 语句用来告诉Python跳过当前循环的剩余语句,然后继续进行下一轮循环。
                                #continue 语句跳出本次循环,而break跳出整个循环。
                        lines = line.strip("\n").split("\t")
                        chr = lines[0] #染色体
                        type = lines[2] #如:gene,exon,...
                        start = int(lines[3])
                        end = int(lines[4])
                        orientation = lines[6]
                        attr = lines[8]

                        #m = re.search(pattern, string)  搜索整个字符串,直到发现符合的子字符串。
                        #import re
                        #m = re.search('[0-9]','abcd4ef')
                        #print(m.group(0))
                        #re.search()接收两个参数,第一个'[0-9]'表示从字符串想要找的是从0到9的一个数字字符;
                        #re.search()如果从第二个参数找到符合要求的子字符串,就返回一个对象m,
                        #通过m.group()的方法查看搜索到的结果。如果没有找到符合要求的字符,re.search()会返回None。
                        
                        if not chr in list_chr:
                                list_chr.append(chr) 

                        if type == "gene":
                                gene = Gene()
                                id = re.search(r'ID=([^;]+);?',attr).group(1)
                                #[url=http://www.runoob.com/python/python-reg-expressions.html]http://www.runoob.com/python/python-reg-expressions.html[/url]
                                ##group(0)是整个正则表达的搜索结果,group(1)是第一个群……
                                #"+"      匹配1个或多个的表达式
                        #"?"      匹配0个或1个由前面的正则表达式定义的片段
                        #[^...]          不在[]中的字符; [^abc] 匹配除了a,b,c之外的字符。
                        #(re)          匹配括号内的表达式,也表示一个组
                                gene.chr = chr
                                gene.start = start
                                gene.end = end
                                gene.orientation = orientation
                                list_gene(id) = gene
                                #print(id)
                        elif type == "mRNA":
                                transcript = Transcript()
                                id = re.search(r'ID=([^;]+);?',attr).group(1)
                                parent =re.search(r'Parent=([^;]+);?',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'Parent=([^;]+);?',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):
        """
        chr_gene
        :param list_gene:
        :return:
        """

        print("chr_gene")
        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):
        """
        gene_len
        :param list_gene:
        :return:
        """

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

def gene_transcript(list_transcript):
        """
        gene_transcript
        :param list_transcript:
        :return:
        """

        print("gene_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 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):
        """
        transcript_exon
        :param list_exon:
        :return:
        """

        print("transcript_exon") 
        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):
        """
        exon_pos
        :param list_exon:
        :return:
        """

        print("exon_pos")
        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")


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



3. 第四题和第三题均使用python完成,并且,根据教学视频和不同的实验材料,对代码进行适当的思考与调整,从而能更好的理解整个程序,同时也觉得自己python的能力得到了大大的提高,感谢东哥!!!

回复 支持 反对

使用道具 举报

0

主题

5

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-2-16 22:12:43 | 显示全部楼层
Perl-182-
首先,下载人类基因组Gene set文件,产看文件格式
可以看到,数据共分为9列,每一列的意义可以根据说明文件进行了解:


#查看文件是用到zcat命令:用于不真正解压压缩文件就可以显示压缩文件中的内容。同时我们只看前几行就行,所以命令为:
zcatHomo_sapiens.GRCh38.87.chr.gtf.gz | head -10 #查看前10行,如果查看1000-1005行,可以写成:zcat Homo_sapiens.GRCh38.87.chr.gtf.gz| Perl –e ‘{while(<>){print if (1000..1005)}}’
##############################################################################
结题思路:如果统计每条染色体上基因的数目,那么只需要将第三列为“gene”的行取出,然后统计第一列(染色体)对应多少个基因即可。
[Perl] 纯文本查看 复制代码
open (DATA, 'gzip -dc Homo_sapiens.GRCh38.87.chr.gtf.gz|') or die ("can not open Homo_sapiens.GRCh38.87.chr.gtf.gz");

while(<DATA>){ #读取数据,在读取数据之前进行解压
	next if (/^#!/); #跳过开头为#!的行
	@lines = split /\t/, $_; #对每一行进行split,然后就结果存入数组@lines
		if ($lines[2] eq "gene"){ #统计数组第三位为“gene”的行
			$chr = $lines[0]; #将染色体号赋值给变量$chr
			$chr_gene_number{$chr} +=1; #哈希,这是本题的关键!!!键值为染色体号,键值相同的,就对对应的染色体基因数目加1
		}
	}
close DATA;

foreach $chr(sort keys %chr_gene_number){ #对哈希键进行排序
print "$chr\t$chr_gene_number{$chr}\n"; #输出键(染色体号)及其对应的基因数
}

运行结果:

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

4

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2017-2-19 01:13:52 | 显示全部楼层
我是做牛的,下载了牛Btau_5.0.1的基因组,我的存在一个特殊情况就是gene_biotype可能在行末出现,而我下载的基因组行末没有分号(;),所以我首先在行末加了分号s/$/;/之后进行的统计。

#!/usr/bin/perl
use strict;
my @gene;
my @transcript;
while (<>){
        chomp;
        next if (/^#/);
        s/$/;/;
        my @temp=split /\s/,$_;
        if ($temp[2] eq "gene" ){
                push @gene,$temp[0];
                /gene_biotype=(.*?);/;
                push @transcript,$1;
        }
}
#print "$gene[1]\n";
#print "$transcript[1]\n";
my $gene;
my $transcript;
my %count_gene;
my %count_transcript;
foreach $gene (@gene){
        $count_gene{$gene}++;
}
foreach $transcript (@transcript){
        $count_transcript{$transcript}++
}
foreach $gene (sort keys %count_gene){
        print "$gene\t$count_gene{$gene}\n";
}
foreach $transcript (sort keys %count_transcript){
        print "$transcript\t$count_transcript{$transcript}\n";
}

结果如下;
NC_007299.6     1266
NC_007300.6     1237
NC_007301.6     1833
NC_007302.6     1136
NC_007303.6     1709
NC_007304.6     948
NC_007305.6     1806
NC_007306.6     1102
NC_007307.6     817
NC_007308.6     1450
NC_007309.6     1310
NC_007310.6     662
NC_007311.6     1104
NC_007312.6     734
NC_007313.6     1485
NC_007314.5     947
NC_007315.6     862
NC_007316.6     1657
NC_007317.6     1569
NC_007318.6     512
NC_007319.6     854
NC_007320.6     764
NC_007324.6     1136
NC_007325.6     473
NC_007326.6     972
NC_007327.6     560
NC_007328.5     379
NC_007329.6     442
NC_007330.6     912
NC_007331.5     1647
NC_016145.1     297
NW_014647391.1  1
NW_014647392.1  1
NW_014647411.1  1
NW_014647427.1  1
NW_014647882.1  1
NW_014648117.1  1
NW_014648401.1  1
NW_014648481.1  4
NW_014648483.1  1
NW_014648484.1  4
NW_014648489.1  4
NW_014648492.1  1
NW_014648493.1  2
NW_014648494.1  3
NW_014648496.1  1
NW_014648497.1  1
NW_014648504.1  1
NW_014648506.1  1
NW_014648689.1  1
NW_014648735.1  1
NW_014648766.1  1
NW_014649280.1  1
NW_014649281.1  1
NW_014649366.1  1
NW_014649370.1  1
NW_014649373.1  1
NW_014649411.1  1
NW_014649455.1  2
NW_014649485.1  1
NW_014649525.1  1
NW_014649527.1  1
NW_014649542.1  1
NW_014649547.1  1
NW_014649548.1  2
NW_014649556.1  1
NW_014649563.1  2
NW_014649592.1  1
NW_014649595.1  2
NW_014649611.1  1
NW_014649624.1  1
NW_014649625.1  2
NW_014649681.1  1
NW_014649709.1  1
NW_014649719.1  3
NW_014649726.1  1
NW_014649763.1  1
NW_014649784.1  1
NW_014649786.1  1
NW_014649787.1  1
NW_014649794.1  1
NW_014649799.1  1
NW_014649801.1  1
NW_014649815.1  1
NW_014649816.1  1
NW_014649844.1  1
NW_014649851.1  1
NW_014649853.1  1
NW_014649854.1  1
NW_014649859.1  1
NW_014649924.1  1
NW_014649943.1  1
NW_014649956.1  2
NW_014649973.1  3
NW_014649982.1  1
NW_014649985.1  1
NW_014649988.1  1
NW_014649990.1  1
NW_014649997.1  1
NW_014649998.1  1
NW_014650017.1  1
NW_014650023.1  1
NW_014650026.1  1
NW_014650055.1  1
NW_014650070.1  1
NW_014650095.1  2
NW_014650098.1  1
NW_014650101.1  1
NW_014650156.1  2
NW_014650157.1  1
NW_014650158.1  1
NW_014650159.1  1
NW_014650162.1  1
NW_014650167.1  2
NW_014650174.1  1
NW_014650176.1  1
NW_014650180.1  1
NW_014650267.1  1
NW_014650282.1  1
NW_014650376.1  1
NW_014650421.1  1
C_region        20
RNase_MRP_RNA   1
SRP_RNA 1
V_segment       180
lncRNA  4647
miRNA   806
misc_RNA        162
other   99
protein_coding  21427
pseudogene      3768
rRNA    2
tRNA    1583
telomerase_RNA  1

当然我没来得及替换染色体的编号。
回复 支持 反对

使用道具 举报

0

主题

6

帖子

89

积分

注册会员

Rank: 2

积分
89
发表于 2017-2-19 19:45:30 | 显示全部楼层
本帖最后由 banapi 于 2017-2-19 21:42 编辑

158gtf文件
基本结构:染色体编号,数据库来源,内容,起始和终止位置。。。
最开始是#开头的信息


shell
统计整理表格化文件,感觉shell好方便管道符“|”将两个命令隔开,管道符左边命令的输出就会作为管道符右边命令的输入
less -S 文件名 使表格规范输出,-S如果有解压可以解压
cat不可直接展示压缩文件,zcat 文件名 |less -S
grep -v "^#" | less -S 不要开头是#号的,不加v则是开头是#号的,grep是一个强大的文本搜索工具,有很多参数,功能丰富

gene统计
[Shell] 纯文本查看 复制代码
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |grep -v "^#"|awk '$3~/gene/{print $1,"\t",$3}'|sort|uniq -c|less -S
部分结果
   2204 10       gene
   3235 11       gene
   2940 12       gene
   1304 13       gene
   2224 14       gene
   2152 15       gene
   2511 16       gene
   2995 17       gene
   1170 18       gene
transcript统计
[Shell] 纯文本查看 复制代码
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |grep -v "^#"|awk '$3~/transcript/{print $0}'|awk '$0~/protein_coding/{print $1,"\t",$3}'|sort|uniq -c|less -S

部分结果
4421 10       transcript
   9633 11       transcript
   4421 10       transcript
   9633 11       transcript
   8925 12       transcript
   1868 13       transcript
   5446 14       transcript
   5258 15       transcript
   7707 16       transcript
  10020 17       transcript
   2453 18       transcript
  10680 19       transcript






回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2017-3-15 10:37:36 | 显示全部楼层
bioconda暂时没有win版,这道题,我偷个懒,学习所得是https://genome.ucsc.edu/FAQ/FAQformat.html#format4这个链接。
回复 支持 反对

使用道具 举报

0

主题

3

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-4-13 21:50:33 | 显示全部楼层
本帖最后由 幻灵风舞 于 2017-4-13 21:52 编辑

179-python-沈宇东
学习了大神的python编程,自己也手打了一遍,运行成功,运行的数据是tmp文件。通过练习理顺了编程思路,感觉这个程序已经可以很好的拆分数据用于后续分析。

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

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):
    list_chr = [] #染色体的列表
    dict_gene = {} #{基因id:基因实例(含各种信息:染色体信息、基因id、基因起止位点信息)}
    dict_transcript = {} #{转录本id:转录本实例(含各种信息:染色体信息、基因id、转录本id、转录本起止位点信息)}
    list_exon = [] #外显子实例列表(含各种信息:染色体信息、转录本id、外显子起止位点信息)
    with open('tmp','r') 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 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) # .group(1)?
                gene.chr = chr
                gene.start = start
                gene.end = end
                gene.id = id
                gene.orientation = orientation
                dict_gene[gene.id] = gene

            elif type == "transcript":
                transcript = Transcript()
                id = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
                #正则表达式取括号里的东西,后面;?表示匹配一个;
                parent = re.search(r'gene_id "([^;]+)";?',attr).group(1)
                if parent not in dict_gene:
                    #zidian?
                    continue
                transcript.chr = chr
                transcript.parent = parent
                transcript.id = id
                transcript.start = start
                transcript.end = end
                dict_transcript[transcript.id] = transcript

            elif type == "exon":
                exon = Exon()
                parent = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
                if parent not in dict_transcript:
                    continue
                exon.chr = chr
                exon.parent = parent
                exon.start = start
                exon.end = end
                list_exon.append(exon)
 
    chr_gene(dict_gene)
    gene_len(dict_gene)
    gene_transcript(dict_transcript)
    transcript_exon(list_exon)
    exon_pos(list_exon)
    
def chr_gene(dict_gene):
    print("染色体上基因数量分布")  
    #染色体上基因数量分布的计算并写入文件,使用了两个循环,是因为一个是对所有的基因进行数量统计,需要对整个字典进行统计才行。而不能统计一次写入文件一次。         
    count_gene = {}
    for info in dict_gene.values():  
        #dict_gene 在主函数中 dict_gene.values()为基因id,对基因id进行遍历
        chr = info.chr #基因id的染色体号
        if chr not in count_gene: 
            #看染色体是否在count_gene字典里,如果不在字典里,则在定义count_gene字典中定义当前染色体为键,并将其对应的值加1
            count_gene[chr] =1
        else:
            count_gene[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_transcript(dict_transcript):
    print("基因的转录本数量分布")
    count_transcript = {}
    for info in dict_transcript.values(): 
        #dict_gene 在主函数中
        gene_id = info.parent
        if gene_id not 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):
    print("转录本外显子数量统计")
    count_exon = {}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id not in count_exon:
            count_exon[transcript_id] = 1
        else:
            count_exon[transcript_id] += 1
            
    with open("transcript_exon",'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 gene_len(dict_gene):
    print("基因长度分布情况")
    with open("gene_len.txt",'w') as fp_out:
        for gene_id,info in dict_gene.items():
            len = info.end-info.start+1
            print("\t".join([gene_id,str(len)])+"\n")
            fp_out.write("\t".join([gene_id,str(len)])+"\n")
            
def exon_pos(list_exon):
    print("外显子坐标统计")
    count_exon_pos={}
    for exon in list_exon:
        transcript_id = exon.parent
        if transcript_id not in count_exon_pos:
            count_exon_pos[transcript_id] = "%s-%s"%(str(exon.start),str(exon.end))
        else:
            count_exon_pos[transcript_id] += ".%s-%s"%(str(exon.start),str(exon.end))
            
    with open("exon_pos.txt",'w') as fp_out:
        for transcript_id,exon_pos in count_exon_pos.items():
            print("\t".join([transcript_id, exon_pos])+"\n")
            fp_out.write("\t".join([transcript_id, exon_pos])+"\n")
                
if __name__ == "__main__":
    main(args)

回复 支持 反对

使用道具 举报

0

主题

4

帖子

77

积分

注册会员

Rank: 2

积分
77
发表于 2017-7-26 13:46:44 | 显示全部楼层
lilith098 发表于 2017-1-17 13:24
137+python
题目:简单输出染色体中gene、transcript、exon、cds的数量。
状态:看了视频,根据视频修改了 ...

你好,我本来是用perl的,最近在学python,感觉你这个代码写的挺简洁的,于是照着写了一个,但是发现你的脚本里有点逻辑问题,第23行不应该用循环去叠加吧?此处应该用if 判断,去赋值和叠加:
for feature_id in list_feature:
                       chr_gene[chr_id][feature_id] += 1
这样出来的结果就对了。

回复 支持 反对

使用道具 举报

1

主题

10

帖子

80

积分

注册会员

Rank: 2

积分
80
发表于 2018-5-1 19:38:25 | 显示全部楼层
13235695523 发表于 2017-1-12 22:39
上面的代码是python统计每个染色体中gene,exon,transcirpt的个数,感谢东哥的视频。 ...

请问在哪可以看到东哥的视频?
回复 支持 反对

使用道具 举报

0

主题

4

帖子

91

积分

注册会员

Rank: 2

积分
91
发表于 2018-5-9 22:58:01 | 显示全部楼层
033-Perl-R-Pyth 发表于 2017-1-26 17:46
看视频前
思路:
用hash区别不同的chr,用split函数选取目标列,遍历hash输出[mw_shl_code=perl,true]

这个视频有回放吗
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 15:45 , Processed in 0.059627 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.