搜索
楼主: Jimmy

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

  [复制链接]

0

主题

16

帖子

149

积分

注册会员

Rank: 2

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

还是炒的别人的。。。各位牛人的代码我这种小透明理解起来真费劲啊,
[mw_shl_code=python,true]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] = 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] += 1
               
time_end = time.clock()[/mw_shl_code]

回复 支持 反对

使用道具 举报

0

主题

4

帖子

50

积分

注册会员

Rank: 2

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

103-python-如梦

[mw_shl_code=python,true]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)[/mw_shl_code]
回复 支持 反对

使用道具 举报

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
自己编程:照着直播里面的稍微修改了一下 发现竟然可以运行!看来我抄的还是可以,至少没有抄错了!哈哈!!
[mw_shl_code=applescript,true]#!/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";
}[/mw_shl_code]
各个染色体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
暂时只把这个做通了,剩下部分等学习明白了,在写!
吃完午饭,回来研究了一会代码,发现昨天写的代码,把{}放错地方了,真是尴尬!下面是看了直播的视频后,自己参考着写的!
[mw_shl_code=applescript,true]#!/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";
}[/mw_shl_code]
##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的模块。

我的代码:
[mw_shl_code=python,true]# 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: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")
[/mw_shl_code]

代码使用: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的思想再来搞定。




回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-19 20:09:35 | 显示全部楼层
023R
genomicfeatures这个包太过强大,我需要研究研究。
视频中的操作我基本看懂了,还需练习一下。视频中代码贴出来,有问题请大家指出来
[mw_shl_code=applescript,true]#这是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()[/mw_shl_code]
回复 支持 反对

使用道具 举报

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-涤生
[mw_shl_code=shell,true]zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{print if $F[2] eq "gene" }' |cut -f 1 |sort |uniq –c[/mw_shl_code]
[mw_shl_code=shell,true]zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}'  |sort |uniq –c[/mw_shl_code]
[mw_shl_code=perl,true]#!/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";
}[/mw_shl_code]

高级部分待更新

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

使用道具 举报

0

主题

8

帖子

179

积分

注册会员

Rank: 2

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

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

贴上已完成部分的代码:
[mw_shl_code=python,true]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:


[/mw_shl_code]

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

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

[mw_shl_code=bash,true]('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)[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

4

帖子

51

积分

注册会员

Rank: 2

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

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


每条染色体基因数:[mw_shl_code=perl,true]#!/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";
}[/mw_shl_code]
基因类型:[mw_shl_code=applescript,true]#!/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";
}[/mw_shl_code]

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

[mw_shl_code=applescript,true]#!/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";
}[/mw_shl_code]




回复 支持 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, 2020-4-6 16:19 , Processed in 0.028827 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.