搜索
楼主: Jimmy

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

  [复制链接]

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-1-17 18:10:14 | 显示全部楼层
本帖最后由 xxnku 于 2017-2-19 20:21 编辑

还是炒的别人的。。。各位牛人的代码我这种小透明理解起来真费劲啊,
[Python] 纯文本查看 复制代码
import time
from collections import OrderedDict

time_start = time.clock()

dict = OrderedDict()
list_chr = [] #用不到吧?
list_feature = ['gene','transcript','exon','CDS']

with open(r'D:\biotrainee\third\Homo_sapiens.GRCh38.87.chr.gtf','r') as f:
    for line in f:
        if line.startswith('#'):
            continue
        line = line.rstrip().split('\t')
        
        chr_id = line[0]
        if chr_id not in dict.keys():
            #list_chr.append(chr_id)  #用不到
            dict[chr_id] = {}        #字典里面又搞一个字典?
            #dict[chr_id]['gene'] = 0
            #dict[chr_id]['transcript'] = 0
            #dict[chr_id]['exon'] = 0
            #dict[chr_id]['CDS'] = 0
            for i in list_feature:   #用一个for循环代替上面的4个赋值。
                dict[chr_id][i] = 0
        
        #if "gene" == line[2]:
            #dict[chr_id]['gene'] += 1
        #if "transcript" == line[2]:
            #dict[chr_id]['transcript'] += 1
        #if "exon" == line[2]:
            #dict[chr_id]['exon'] += 1
        #if "CDS" == line[2]:
            #dict[chr_id]['exon'] += 1
        else:
            for i in list_feature:  # 用一个for循环代替上面的4个递增
                if i == line[2]:
                    dict[chr_id][i] += 1
                
time_end = time.clock()


回复 支持 反对

使用道具 举报

0

主题

4

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-1-17 18:56:48 | 显示全部楼层
本帖最后由 huamixinhua 于 2017-2-5 17:17 编辑

103-python-如梦

[Python] 纯文本查看 复制代码
def tj(lj):

    zd={}
    f= open(lj, 'r')
    for i in f:
        if not i.startswith('#'):
            lb=i.split('\t')
            sj=(lb[0],lb[2])
            if sj in zd:
                zd.update({sj:zd.get(sj)+1})
            else:
                zd.setdefault(sj,1)
    f.close

    zd2={}
    for i in zd.items():
        if str(i[0][0]).isdigit():
            sj='{:0>2}'.format(str(i[0][0]))
        else:
            sj=str(i[0][0])

        if sj in zd2:
            zd2.update({sj:zd2.get(sj)+";"+i[0][1]+":"+str(i[1])})
        else:
            zd2.setdefault(sj,i[0][1]+":"+str(i[1]))
    for i in sorted(zd2.items()):
        print(i)
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-18 10:29:44 | 显示全部楼层
本帖最后由 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

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-1-18 15:43:16 | 显示全部楼层
本帖最后由 旭日早升 于 2017-2-3 19:55 编辑

201-python-无名
本周的作业是探索基因组注释gtf文件,首先需要对gtf文件有一定了解,我最早是从gencode开始了解的,算是有一定基础了。
我的代码思路:我完成的是基本的作业,群主给了3个例子,即统计每条染色体gene的书目,统计每条染色体coding gene的书目,以及每种gene类型的书目。代码核心是初始化3个字典分别将用于存储3个对应的结果,依次读入一行并且判断,如果不是注释行,则对该行分隔,首先判断染色体并且给每个字典里对应染色体赋值,随后判断是否为gene,如果是,则gene字典对应染色体的基因数加一,再判断gene类型,如果类型为protein_coding,则coding字典对应的染色体基因数加一,统计gene类型的字典则是判断gene_biotype并且累加在字典中。
其他注意:除了核心的代码,我还跟大神学会了使用sys和getopt实现了接口的操作。另外gtf下载下来是gz压缩文件,只有43M,解压缩后1.3G,为了使代码直接使用gz文件,使用了gzip的模块。

我的代码:
[Python] 纯文本查看 复制代码
# Description: Cainiao3_v1
# Author: WKL
# Data: 2017-01-18
# 

import re
import os
from collections import OrderedDict
from operator import itemgetter

import sys
import getopt
import gzip

## change the path 
path = "./"
os.chdir(path)

## usage
def usage():
        print("")
        print("Explore genome annotation gtf file.")
        print("uage: python %s -option <argument>" %sys.argv[0])
        print(" -h/--help ")
        print(" -i <STRING> GTF input file.")
        print(" -o <STRING> Output file.")

## deal with options
try: 
        opts, args = getopt.getopt(sys.argv[1:], "i:o:h", ["help"])
except getopt.GetoptError:
        print("ERROR: Get option error.")
        usage()
        sys.exit(2)

for opt, val in opts:
        if opt in ("-h", "--help"):
                usage()
                sys.exit(1)
        else:
                if opt in ("-i"):
                        infile = val 
                if opt in ("-o"):
                        outfile = val

try:
        infile, outfile
except:
        print("ERROR: Missing option.")
        usage()
        sys.exit(2)

CHR = OrderedDict()
CHR_coding = OrderedDict()
CHR_type = OrderedDict()

## open input file 
if re.search(r"\.gz$", infile):
        with gzip.open(infile, "rt") as f:
                for line in f:
                        if not line.startswith("#"):
                                line = line.rstrip().split("\t")
                                anno = line[-1].split(" ")
                                Chr = line[0]
                                Type = line[2]
                                gene_type = anno[9]
                                if Chr not in CHR:
                                        CHR[Chr] = 0
                                        CHR_coding[Chr] = 0
                                if Type == "gene":
                                        CHR[Chr] += 1
                                        if gene_type not in CHR_type:
                                                CHR_type[gene_type] = 1 
                                        else:
                                                CHR_type[gene_type] += 1
                                        if gene_type == "\"protein_coding\";":
                                                CHR_coding[Chr] += 1
else:
        with open(infile, "rt") as f:
                for line in f:
                        if not line.startswith("#"):
                                line = line.rstrip().split("\t")
                                anno = line[-1].split(" ")
                                Chr = line[0]
                                Type = line[2]
                                gene_type = anno[9]
                                if Chr not in CHR:
                                        CHR[Chr] = 0 
                                        CHR_coding[Chr] = 0 
                                if Type == "gene":
                                        CHR[Chr] += 1
                                        if gene_type not in CHR_type:
                                                CHR_type[gene_type] = 1 
                                        else:
                                                CHR_type[gene_type] += 1
                                        if gene_type == "\"protein_coding\";":
                                                CHR_coding[Chr] += 1

with open(outfile,"wt") as wf:
        wf.write("CHR_gene\n")
        for key, value in CHR.items():
                wf.write(str(key) + "\t" + str(value) + "\n")
        wf.write("\nCHR_coding\n")
        for key, value in CHR_coding.items():
                wf.write(str(key) + "\t" + str(value) + "\n")
        wf.write("\nCHR_type\n")
        for key, value in CHR_type.items():
                wf.write(str(key) + "\t" + str(value) + "\n")


代码使用:python Cainiao_3.py -i Homo_sapiens.GRCh38.87.chr.gtf.gz -o Cainiao_3.output
运行结果:
CHR_gene
1       5194
2       3971
3       3010
。。。
CHR_coding
1       2052
2       1255
3       1076
。。。
CHR_type
"transcribed_unprocessed_pseudogene";   735
"unprocessed_pseudogene";       2661
"miRNA";        1567
。。。
接下来学习老师的视频,再来总结。还是要像群主说的,多找出有意思的问题,多多实践


年后回来又看了东神的视频,数据的结构化确实给了我很大的启发,这种层级衍生的代码的可移植性和操作性确实强大啊,最后东老师留的作业我就有点懵了,一时不知道如何下手啊,先放两天学习下json的思想再来搞定。




回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-19 20:09:35 | 显示全部楼层
023R
genomicfeatures这个包太过强大,我需要研究研究。
视频中的操作我基本看懂了,还需练习一下。视频中代码贴出来,有问题请大家指出来
[AppleScript] 纯文本查看 复制代码
#这是Hg19版本的txdb文件
source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
#我们希望使用Hg38版本,先转换一下
library(GenomicFeatures)
setwd("C:\\Users\\qin\\Desktop")
GRCh38_txdb<-makeTxDbFromGFF("Homo_sapiens.GRCh38.87.chr.gtf.gz",organism = "Homo sapiens")#
#查看变量和染色体
columns(GRCh38_txdb)
seqlevels(GRCh38_txdb)
saveDb(GRCh38_txdb,file = "GRCh38.sqlite")#保存为本地
GRCh38_txdb<-loadDb("GRCh38.sqlite")#载入本地文件

#查看各条染色体的基因分布
grch38_gene<-genes(GRCh38_txdb,columns="gene_id")
grch38_gene_df<-as.data.frame(grch38_gene)
head(grch38_gene_df)
nrow(grch38_gene_df)
#[1] 57992
library(ggplot2)
#geom_bar是对数目的作图,查看帮助指导如果对变量作图用geom_col
ggplot(grch38_gene_df,aes(x=factor(seqnames)))+geom_bar()

#同理查看各转录本
grch38_tx<-transcripts(GRCh38_txdb,columns="tx_name")
grch38_tx_df<-as.data.frame(grch38_tx)
nrow(grch38_tx_df)
#[1] 197935
ggplot(grch38_tx_df,aes(x=factor(seqnames)))+geom_bar()#得到转录本按染色体的分布

#在查看各个基因所含转录本的平均水平
GRCh38_tx_by_genne<-transcriptsBy(GRCh38_txdb,by="gene")
GRCh38_tx_by_genne<- as.data.frame(GRCh38_tx_by_genne)
nrow(as.data.frame(GRCh38_tx_by_genne))
#[1] 197935
library(dplyr)
GRCh38_tx_by_genne<-GRCh38_tx_by_genne%>%group_by(group_name)%>%mutate(count=n())
head(GRCh38_tx_by_genne)
sort(unique(GRCh38_tx_by_genne$count))
summary(GRCh38_tx_by_genne$count)
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1.00    3.00    9.00   11.04   15.00  170.00 
df<-unique(GRCh38_tx_by_genne[,c("group_name","count")])
ggplot(df,aes(x=factor(count)))+geom_bar()

#考察转录本所含外显子数目问题
GRCh38_tx_by_exon<-transcriptsBy(GRCh38_txdb,by="exon")
GRCh38_tx_by_exon<- as.data.frame(GRCh38_tx_by_genne)
nrow(as.data.frame(GRCh38_tx_by_exon))
#[1] 1181908
GRCh38_tx_by_exon<-as.data.frame(GRCh38_tx_by_exon)
GRCh38_tx_by_exon<-GRCh38_tx_by_exon%>%group_by(tx_name)%>%mutate(count=n())
summary(GRCh38_tx_by_exon$count)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1.00    5.00    9.00   13.68   17.00  363.00 
df_exon<-unique(GRCh38_tx_by_exon[,c("tx_name","count")])
ggplot(df_exon,aes(x=factor(count)))+geom_bar()
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-19 23:41:47 | 显示全部楼层
{summary(GRCh38_tx_by_exon$count)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
1.00    5.00    9.00   13.68   17.00  363.00}这一步好像不对,中位数是9简直了
得到数据框后再summary就没问题
summary(df_exon)
group_name            count        
Length:57992       Min.   :  1.000  
Class :character   1st Qu.:  1.000  
Mode  :character   Median :  1.000  
                    Mean   :  3.413  
                    3rd Qu.:  4.000  
                    Max.   :170.000  
平均每个基因有3个左右的转录本,符合我的印象。
回复 支持 反对

使用道具 举报

1

主题

15

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-21 18:42:20 | 显示全部楼层
本帖最后由 disheng 于 2017-1-22 19:13 编辑

043-Perl+R+shell-涤生
[Shell] 纯文本查看 复制代码
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{print if $F[2] eq "gene" }' |cut -f 1 |sort |uniq –c

[Shell] 纯文本查看 复制代码
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}'  |sort |uniq –c

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

my %gene;
while(<>){
        s/^#.*\n$//g;
        @lines = split /\t/;
        if($lines[3] eq "gene"){
                $gene{$lines[0]}++;
        }
}
foreach $chr( sort keys %gene){
        print "$chr\t$gene{$chr}\n";
}


高级部分待更新

第三国.png
韩国返回法国.png
回复 支持 反对

使用道具 举报

0

主题

8

帖子

179

积分

注册会员

Rank: 2

积分
179
发表于 2017-1-22 21:53:51 | 显示全部楼层
028-python-Ryu

因为这周连续上班,今天周日还要回公司上班,打算先凭自己能力先完成作业,再看老师的视频改进自己的程序,因此作业还没有完成,对老师们说声抱歉。但是基础作业的大概框架已经做完,有时间会尽快补上。

贴上已完成部分的代码:
[Python] 纯文本查看 复制代码
import gzip
import re

gtf_file = "Homo_sapiens.GRCh38.87.chr.gtf.gz"

with gzip.open(gtf_file, 'r') as human_gtf:

    gene_dict={} #NOTE: {chr:[{gene1:coordinate1},{gene2:coordinate2}]}
    transcript_dict={} #NOTE: gene:transcript:coordinate
    exon_dict={} #NOTE: transcript:exon:coordinate
    cds_dict={} #NOTE: tarnascript:cds:coordinate
    for line in human_gtf:
        if "#!" not in line:
            content = line.strip().split("\t")
            # print(content)
            
            chr=content[0]
            coordinate=chr+":"+content[3]+"-"+content[4]

            gene_id = re.findall(r'gene_id "(.*?)";',content[8])[0] #NOTE: gene id must be found,so [0]
            transcript_id = re.findall(r'transcript_id "(.*?)";',content[8])
            exon_number = re.findall(r'exon_number "(.*?)";',content[8])
            cds_number = re.findall(r'exon_number "(.*?)";',content[8])


            # TITLE: generate data sets
            ## TITLE: analysis gene
            if content[2] == "gene":
                # print(gene_id)
                if chr in gene_dict:
                    gene_dict[chr].append({gene_id:coordinate})
                else:
                    gene_dict[chr]=[{gene_id:coordinate}]
                # print(gene_dict)

            ## TITLE: analysis transcript
            elif content[2] == "transcript":
                if gene_id in transcript_dict:
                    transcript_dict[gene_id].append({transcript_id[0]:coordinate})
                else:
                    transcript_dict[gene_id]=[{transcript_id[0]:coordinate}]
            
            ## TITLE: analysis exon_dict
            elif content[2] == "exon":
                if transcript_id[0] in exon_dict:
                    exon_dict[transcript_id[0]].append({exon_number[0]:coordinate})
                else:
                    exon_dict[transcript_id[0]]=[{exon_number[0]:coordinate}]
            
            elif content[2] == "CDS":
                if transcript_id[0] in cds_dict:
                    cds_dict[transcript_id[0]].append({cds_number[0]:coordinate})
                else:
                    cds_dict[transcript_id[0]]=[{cds_number[0]:coordinate}]

    for k in gene_dict:
        print(k,len(gene_dict[k]))
    # for k,v in transcript_dict:




目前代码完成情况是这样的,
1.已经将Ensembl的Human的GTF注释文件解析,并且分为基因、转录本、外显子、CDS四个字典,存储了它们的id与范围
2.目前已经完成统计每条染色体上的基因个数
3.打算下一步统计基因、转录本、外显子、CDS的平均个数以及长度范围分布等,并且使用python绘图
4.目前的代码在Ubuntu上我的机子运行51s左右(win上速度会比ubuntu慢),直接读取gzip文件,因为下一步的多物种的GTF统计如果全部解压一遍实在麻烦而且占空间。

将目前的基因个数的统计结果附上,作业完成后更新此贴

[Bash shell] 纯文本查看 复制代码
('20', 1386)
('21', 835)
('22', 1318)
('1', 5194)
('3', 3010)
('2', 3971)
('5', 2868)
('4', 2505)
('7', 2867)
('6', 2863)
('9', 2242)
('8', 2353)
('Y', 523)
('X', 2359)
('11', 3235)
('10', 2204)
('13', 1304)
('12', 2940)
('15', 2152)
('14', 2224)
('17', 2995)
('16', 2511)
('19', 2926)
('18', 1170)
('MT', 37)
回复 支持 反对

使用道具 举报

0

主题

4

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-1-23 13:30:06 | 显示全部楼层
本帖最后由 zhongzheng 于 2017-1-23 14:16 编辑

212-perl-zhongzheng
对视频中给出的perl代码做了些修改,运行结果一致。


每条染色体基因数:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use warnings;
use strict;

my ($row,@cells,@search_feature,$search_feature,@search_cells,@chr,$count,%count,$chr);
open DATA, "Homo_sapiens.GRCh38.87.chr.gtf";
open OUTPUT, ">gene_number.txt";

while ($row=<DATA>) {
    last unless $row =~ /\S/;
        @cells = split (/\t/, $row);
        if ($cells[2] eq "gene"){
                $count{$cells[0]}++;
                }
}

foreach $chr(sort keys %count){
    print OUTPUT "$chr\t=>\t$count{$chr}\n";
}

基因类型:
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/perl
use warnings;
use strict;

my ($row,@cells,$gene_biotype,@gene_biotype,$count,%count);
open DATA, "Homo_sapiens.GRCh38.87.chr.gtf";
open OUTPUT, ">gene_biotype.txt";

while ($row=<DATA>) {
last unless $row =~ /\S/;
@cells = split (/\t/, $row);
if ($cells[2] eq "gene") {
        if ($row =~ /gene_biotype "(.*?)";/){
            $count{$1}++;
        } 
    }
}

foreach $gene_biotype(sort keys %count){
print OUTPUT "$gene_biotype\t=>\t$count{$gene_biotype}\n";
}


每条染色体蛋白编码基因数:

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

my ($row,@cells,$gene_biotype,@gene_biotype,$count,%count);
open DATA, "Homo_sapiens.GRCh38.87.chr.gtf";
open OUTPUT, ">protein_coding.txt";

while ($row=<DATA>) {
    last unless $row =~ /\S/;
        @cells = split (/\t/, $row);
        if ($cells[2] eq "gene") {
            if ($row =~ /gene_biotype "protein_coding";/){
                        $count{$cells[0]}++;
                } 
        }
}

foreach $row(sort keys %count){
    print OUTPUT "$row\t=>\t$count{$row}\n";
}





回复 支持 0 反对 2

使用道具 举报

0

主题

4

帖子

49

积分

新手上路

Rank: 1

积分
49
发表于 2017-1-24 10:01:24 | 显示全部楼层
用R
source("http://bioconductor.org/biocLite.R")
biocLite("GenomicFeatures")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
browseVignettes(package = "TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb<- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb
library("GenomicFeatures")
GRCh38_txdb<-makeTxDbFromGFF("Homo_sapiens.GRCh38.87.chr.gtf.gz",
                             dataSource = "ENSEMBL",
                             organism = "Homo sapiens")
saveDb(GRCh38_txdb,file="GRCh38.87.sqlite")
GRCh38_txdb<-loadDb(GRCh38.87.sqlite)
columns(GRCh38_txdb)
seqlevels(GRCh38_txdb)

grch38_genes<-genes(GRCh38_txdb,columns="gene_id")
head(as.data.frame(grch38_genes))
grch38_genes_df<-as.data.frame(grch38_genes)
nrow(grch38_genes_df)
library(ggplot2)
ggplot(grch38_genes_df,aes(x=factor(seqnames)))+geom_bar()
#以上为每条染色体基因个数
grch38_tx<-transcripts(GRCh38_txdb,columns=c("tx_name"))
head(as.data.frame(grch38_tx))
grch38_tx_df<-as.data.frame(grch38_tx)
nrow(grch38_tx_df)
ggplot(grch38_tx_df,aes(x=factor(seqnames)))+geom_bar()
table(grch38_genes_df$seqnames)
#以上为每条染色转录本个数

grch38_tx_by_genes<-transcriptsBy(GRCh38_txdb,by=c("gene"))
head(as.data.frame(grch38_tx_by_genes))
nrow(as.data.frame(grch38_tx_by_genes))
grch38_tx_by_genes<-as.data.frame(grch38_tx_by_genes)
head(grch38_tx_by_genes)
library(dplyr)
grch38_tx_by_genes<-grch38_tx_by_genes %>% group_by(group_name) %>% mutate(count=n())
head(grch38_tx_by_genes)
sort(as.numeric(unique(grch38_tx_by_genes$count)))
summary(as.numeric(grch38_tx_by_genes$count))
df<-unique(grch38_tx_by_genes[,c("group_name","count")])
ggplot(df,aes(x=factor(count)))+geom_bar()
head(df)
#转录本比基因

grch38_exons<-transcriptsBy(GRCh38_txdb,by=c("exon"))
head(as.data.frame(grch38_exons),10)
nrow(as.data.frame(grch38_exons))
grch38_exons<-as.data.frame(grch38_exons)
grch38_exons<-grch38_exons %>% group_by(tx_name) %>% mutate(count=n())
head(grch38_exons)
sort(as.numeric(unique(grch38_exons$count)))
summary(as.numeric(grch38_exons$count))
df_exons<-unique(grch38_exons[,c("tx_name","count")])
ggplot(df_exons,aes(x=factor(count)))+geom_bar()
#转录本比exon
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-20 21:07 , Processed in 0.037488 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.