搜索
楼主: Jimmy

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

  [复制链接]

0

主题

4

帖子

97

积分

注册会员

Rank: 2

积分
97
发表于 2017-2-10 19:05:50 | 显示全部楼层
本帖最后由 youzicha 于 2017-2-13 14:37 编辑

#编号099
老师讲的很细致,我把代码整合到一个里面了,如下
[Perl] 纯文本查看 复制代码
#! /usr/bin/perl 
use warnings;
use strict;
my (%count,$line,@line,$item,%count_t,%count_p,$i,$items);
################# select gene and count ###############
open (FILE,"<Homo_sapiens.GRCh38.87.chr.gtf")||die "$!";
open (OUT,">gene_count.txt")||die "$!";
while ($line=<FILE>){
        next if ($line=~/^#/);
        @line=split (/\t/,$line);
        if ($line[2]=~/gene/){
                $count{$line[0]}++;
        }
}

foreach $item( sort keys %count){
        print OUT "$item\t => \t $count{$item}\n";
}
close OUT;
############### select type and count #################
open (TYPE,">type_count.txt")||die "$!";
open (FILE,"<Homo_sapiens.GRCh38.87.chr.gtf")||die "$!";
while ($line=<FILE>){
        next if ($line=~/^#/);
        @line=split (/\t/,$line);
        if ($line[2]=~/gene/){
                if ($line[8]=~/gene_biotype "(.*?)";/){
                        $count_t{$1}++;
                }
        }
}
foreach $i( sort keys %count_t){
        print TYPE "$i\t => \t $count_t{$i}\n";
}
close TYPE;
############### select protein_coding and count #################
open (FILE,"<Homo_sapiens.GRCh38.87.chr.gtf")||die "$!";
open (PRO,">protein_coding_count.txt")||die "$!";
while ($line=<FILE>){
        next if ($line=~/^#/);
        @line=split (/\t/,$line);
        if ($line[2]=~/gene/){
                if ($line[8]=~/gene_biotype "(protein_coding)";/){
                        $count_p{$line[0]}++;;
                }
        }
}
foreach $items( sort keys %count_p){
        print PRO "$items\t => \t $count_p{$items}\n";
}
close PRO;
close FILE;

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

3

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-2-11 00:23:36 | 显示全部楼层

082-python/R-wiwi

082-python/R-wiwi

[AppleScript] 纯文本查看 复制代码
#对gtf文件进行探究
#gtf/gff:[url=http://asia.ensembl.org/info/website/upload/gff.html]http://asia.ensembl.org/info/website/upload/gff.html[/url]



#安装packages:GenomicFeatures\ TxDb.Hsapiens.UCSC.hg19.knownGene

source("https://bioconductor.org/biocLite.R")
biocLite("AnnotationDbi")
source("https://bioconductor.org/biocLite.R")
biocLite("GenomicFeatures")
library("GenomicFeatures")
source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")

#txdb
GRCh38_txdb <- makeTxDbFromGFF("Homo_sapiens.GRCh38.87.chr.gtf.gz",dataSource  ="ENSEMBL",organism="Homo sapiens",format = "gtf")
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)
head(grch38_genes_df)
nrow(grch38_genes_df)
library(ggplot2)
ggplot(grch38_genes_df,aes(x=factor(seqnames)))+geom_bar()



[AppleScript] 纯文本查看 复制代码
#转录本
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()


[AppleScript] 纯文本查看 复制代码
#所有基因平均有多少个转录本
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)
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()


[AppleScript] 纯文本查看 复制代码
#基因上的外显子分布情况?
grch38_exons <- transcriptsBy(GRCh38_txdb,by=c("exon"))
head(as.data.frame((grch38_exons)),10)
grch38_exons <- as.data.frame(grch38_exons)
nrow(grch38_exons)
grch38_exons <- grch38_exons %>% group_by(tx_name) %>% mutate(count=n())
sort(as.numeric(unique(grch38_exons$count)))
summary(as.numeric(grch38_exons$count))
df <- unique(grch38_exons[,c("group_name","count")])
ggplot(df,aes(x=factor(count)))+geom_bar()







回复 支持 反对

使用道具 举报

1

主题

5

帖子

87

积分

注册会员

Rank: 2

积分
87
发表于 2017-2-11 23:11:16 | 显示全部楼层
191-芃


[AppleScript] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my ($hash,($chr, %gene_count,%exon_count);

open(IN,"$ARGV[0]");

while(<IN>){
    chmop;
    if (/^#/){
        next;
      }
      else{
         $chr=(split(/\t/,$_))[0];
         if( (split("\t",$_))[2] eq 'gene'){
             $g_count{$chr}=$g_count{$chr}+1;
             }
                 elsif ( (split("\t",$_))[2] eq 'exon'){
                     $e_count{$chr}=$e_count{$chr}+1;
                }
        }
}
close IN;

foreach (sort keys %g_count){
    print "chr$_\tgene\t$g_count{$_}\n";
    print "chr$_\texon\t$e_count{$_}\n";
}




结果和大家相同,不另付在后面啦
回复 支持 反对

使用道具 举报

0

主题

6

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2017-2-11 23:36:20 | 显示全部楼层
本帖最后由 snjiang12 于 2017-2-12 10:55 编辑

122-R-江
说一下,要注意安装那个GenomicFeatures包的时候,视频里仅仅是在bioconductor网站上复制了代码就完成了,而我个R这根本不行,找了很多帖子,最后发现安装这个GenomicFeatures的时候必须要一些其它的包,这些其它的包也只能通过网上下载后自己解压后复制到library所在的文件中去才行(老是提示版本不支持)。那些很大的包还是自己下载比较好,链接复制到迅雷里去一会就搞完了。
代码搞过来:
GRCh38_tlibrary("GenomicFeatures")#导入GenomicFeatures包
GRCh38_txdb<-makeTxDbFromGFF("Homo_sapiens.GRCh38.87.chr.gtf",dataSourse="ENSEMBL",organism="Homo sapiens")
saveDb(GRCh38_txdb,file="GRCh38.87.sqlite")
GRCh38_txdb<-loadDb("GRCh38.87.sqlite")#在GenomicFeatures当中需要用loadDb函数读取.sqlite文件
columns(GRCh38_txdb)#看下有哪些列
seqlevels(GRCh38_txdb)#染色体名
grch38_genes<-genes(GRCh38_txdb,columns="gene_id")
grch38_genes_df<-as.data.frame(grch38_genes)#变为数据框格式
nrow(grch38_genes_df)
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)
library(ggplot2)
ggplot(grch38_tx_df,aes(x=factor(seqnames)))+geom_bar()#这里用了因子,将变量转换为一系列整数。
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())#%>%相当于增加的意思
sort(as.numeric(unique(grch38_tx_by_genes$count)))
summary(as.numeric(unique(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)
grch38_exons<-as.data.frame(grch38_exons)
grch38_exons<-grch38_exons %>% group_by(tx_name) %>% mutate(count=n())
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()
没啥可说的了,都是理解而已。
回复 支持 反对

使用道具 举报

1

主题

10

帖子

112

积分

注册会员

Rank: 2

积分
112
发表于 2017-2-12 20:33:43 | 显示全部楼层
127
最近,发文章,老板催着补实验补数据,过年也没回家,python没有任何基础,只能零零星星学一点,学多少交多少!
现在卡在正则表达式上了,正在学,后面接着做!现在知道正则表达式为什么要加r,学元字符的含义,\n\r\t的区别,编程目前太弱了~


回复 支持 反对

使用道具 举报

0

主题

2

帖子

111

积分

注册会员

Rank: 2

积分
111
发表于 2017-2-13 10:24:08 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#! usr/bin/perl
use warnings;
use strict;

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

while ( $row = <DATA> ){
  last unless $row =~ /\s/;
  @cells = split /\t/, $row ||next;
 $cells[2]||next;
  if ( $cells[2] eq "gene"){

   if($row =~/gene_biotype(.*?);/){
	
$count{$1}++
	}
  }
}

foreach my $gene_biotype1(sort keys %count){
  printf  "%-40s=>$count{$gene_biotype1}\n",$gene_biotype1;
 }

回复 支持 反对

使用道具 举报

0

主题

2

帖子

111

积分

注册会员

Rank: 2

积分
111
发表于 2017-2-13 10:26:03 | 显示全部楼层
tswcbyolkw 发表于 2017-2-13 10:24
[mw_shl_code=perl,true]#! usr/bin/perl
use warnings;
use strict;

116赞,简单的代码,稍微一点点错误,就不能运行了,检查了半天
回复 支持 反对

使用道具 举报

0

主题

3

帖子

42

积分

新手上路

Rank: 1

积分
42
发表于 2017-2-13 22:07:20 | 显示全部楼层
大致思路是用python中的pandas包构建Dataframe来解决
但是遇到了问题:
1、python3中如何建一个空的,全局的,dataframe
2、pd.merge 几个参数的设置
import os
import pandas as pd
add = os.chdir('C:\\Users\\59779\\bioinfo\\biotree\\4 合并文件\\GSE48213_RAW')
gene_id = []
ENS_index = []
value = []
n = pd.DataFrame()
for file in os.listdir(add):
    if file.endswith('.txt'):
        with open(file) as fh:
            for line in fh:
                if line.startswith('EnsEMBL_Gene_ID'):
                    id = line.split('\t')[1].rstrip()
                    gene_id.append(id)
                    continue
                ENS_index.append(line.split('\t')[0])
                value.append(line.split('\t')[1].rstrip())        
        m = pd.DataFrame(value,index = ENS_index , columns = gene_id)
        n = pd.merge(m,n)
        
print(n)希望有了解的大神能帮助沿着该思路解决一下、
回复 支持 反对

使用道具 举报

0

主题

6

帖子

337

积分

中级会员

Rank: 3Rank: 3

积分
337
发表于 2017-2-14 10:20:13 | 显示全部楼层
本帖最后由 sophieliu 于 2017-2-14 10:24 编辑

035+python+shell
1.每条染色体基因的个数
[Shell] 纯文本查看 复制代码
cat Homo_sapiens.GRCh38.87.chr.gtf | awk '{if($3=="gene"){print $0}}' | awk '{a[$1]++}END{for(i in a){print i,a[i]}}' 


2.每个染色体上protein coding gene统计:
[Shell] 纯文本查看 复制代码
cat Homo_sapiens.GRCh38.87.chr.gtf | awk '$3~/gene/{print $0}' | grep "protein_coding" | awk '{a[$1]++}END{for(i in a){print i,a[i]}}'

3.每个基因的类型统计:
[Shell] 纯文本查看 复制代码
cat Homo_sapiens.GRCh38.87.chr.gtf | awk '$3~/gene/{print $0}' | cut -d ";" -f5 | awk '{a[$2]++}END{for(i in a){print i, a[i]}}'
回复 支持 反对

使用道具 举报

0

主题

5

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2017-2-14 10:29:40 | 显示全部楼层
本帖最后由 死小孩 于 2017-2-14 11:52 编辑

109-R+perl+python-死小孩
[Python] 纯文本查看 复制代码
from collections import OrderedDict

chr_gene = OrderedDict()

#with open ('Homo_GRCh38.87.chr.gtf','r') as gtf:
with open ('Homo_sapiens.GRCh38.87.chr.gtf','r') as gtf:
    for line in gtf:
        num=0
        line=line.rstrip()
        if line.startswith('#'):continue
        chr=line.split('\t')[0]
        if chr not in chr_gene.keys():
            chr_gene[chr]=0
        elif line.split('\t')[2]=='transcript':
            chr_gene[chr]+=1
print(chr_gene)


[Perl] 纯文本查看 复制代码
open PO,"$ARGV[0]" || die "Cannot open hg19_file $ARGV[0] \n";
open OUTFILE,">$ARGV[1]" || die "Cannot open output_file $ARGV[1] \n";

my $line = <PO>;
@line1 = split("\t",$line);
$chr_name = $line1[0];
$gene_number = 1;
while(<PO>){
	$line1[$i] = $_;
	@line_vals1 = split("\t",$line1[$i]);     #按\t对每一行行分隔;
	$chr_id = $line_vals1[0];
	if($chr_name eq $chr_id){
		if ($line_vals1[2] =~ /gene/){
			$gene_number += 1;
		}
		if($line1[$i] =~ /UTR/){
			if($line_vals1[2] =~ /5/){
				$number5 += 1;
			}
			if($line_vals1[2] =~ /3/){
				$number3 += 1;
			}
		}		
		if($line_vals1[2] =~ /exon/){
			$exon_number += 1;
		}	
	} 
	else{
		print OUTFILE "$chr_name\t"."gene\t$gene_number\t"."5UTR\t$number5\t"."exon\t$exon_number\t"."3UTR\t$number3\n";
		$chr_name = $chr_id;
		undef $gene_number;
		undef $number5;
		undef $exon_number;
		undef $number3;
	}
}
print OUTFILE "$chr_id\t"."gene\t$gene_number\t"."5UTR\t$number5\t"."exon\t$exon_number\t"."3UTR\t$number3\n";
		

回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.