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

youzicha

#编号099
老师讲的很细致，我把代码整合到一个里面了，如下[mw_shl_code=perl,true]#! /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;[/mw_shl_code]

寒冰wiwi

082-python/R-wiwi

[mw_shl_code=applescript,true]#对gtf文件进行探究
#gtf／gff：http://asia.ensembl.org/info/website/upload/gff.html



#安装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()[/mw_shl_code]


[mw_shl_code=applescript,true]#转录本
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()
[/mw_shl_code]

[mw_shl_code=applescript,true]#所有基因平均有多少个转录本
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()
[/mw_shl_code]

[mw_shl_code=applescript,true]#基因上的外显子分布情况？
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()
[/mw_shl_code]

吃瓜群众

191-芃


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



结果和大家相同，不另付在后面啦

snjiang12

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()
没啥可说的了，都是理解而已。

鹏洛克贾路

127
最近，发文章，老板催着补实验补数据，过年也没回家，python没有任何基础，只能零零星星学一点，学多少交多少！
现在卡在正则表达式上了，正在学，后面接着做！现在知道正则表达式为什么要加r，学元字符的含义，\n\r\t的区别，编程目前太弱了~

tswcbyolkw

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

tswcbyolkw

tswcbyolkw 发表于 2017-2-13 10:24
[mw_shl_code=perl,true]#! usr/bin/perl
use warnings;
use strict;

116赞,简单的代码，稍微一点点错误，就不能运行了，检查了半天

Irene_xzw

大致思路是用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)希望有了解的大神能帮助沿着该思路解决一下、

sophieliu

035+python+shell
1.每条染色体基因的个数
[mw_shl_code=shell,true]cat Homo_sapiens.GRCh38.87.chr.gtf | awk '{if($3=="gene"){print $0}}' | awk '{a[$1]++}END{for(i in a){print i,a}}' [/mw_shl_code]

2.每个染色体上protein coding gene统计：
[mw_shl_code=shell,true]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}}'[/mw_shl_code]
3.每个基因的类型统计：
[mw_shl_code=shell,true]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}}'[/mw_shl_code]

死小孩

109-R+perl+python-死小孩
[mw_shl_code=python,true]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)[/mw_shl_code]

[mw_shl_code=perl,true]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";
                [/mw_shl_code]