搜索
查看: 555|回复: 2

蛋白编码区域的GC含量会比基因组其它区域的高吗?

[复制链接]

482

主题

848

帖子

2870

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
2870
发表于 2016-10-7 23:06:25 | 显示全部楼层 |阅读模式
可以拿hg19基因组及对应的注释文件来写脚本探究!蛋白编码区域可以从NCBI的CCDS数据库Consensus CDS (CCDS)里面拿到,就是蛋白编码区域,约 32 Mbp, 人基因组 ,hg19,约3G。基因组上面还有很多特殊的区域,比如: centromeres and telomeres 他们的GC含量尤其不一样,有兴趣也可以探究一下。
这样很容易写脚本来统计不同区域的GC含量,可以参考我博客关于CCDS探究的文章:http://www.bio-info-trainee.com/1108.html




上一篇:第二贴:CpG岛
下一篇:关于血液病肿瘤不得不知的两个融合基因(上)
回复

使用道具 举报

5

主题

24

帖子

159

积分

版主

Rank: 7Rank: 7Rank: 7

积分
159
发表于 2016-10-9 17:05:43 | 显示全部楼层
答曰:确实高。下面是我的分析过程:

1、从CCDS官网下载人类CDS序列。 ftp://ftp.ncbi.nih.gov/pub/CCDS/ CDDS的ftp站,点击就能下载

2、编写perl脚本
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w

use strict;
use List::Util qw( sum );
use 5.18.2;

my %data;

while(<>) {
	chomp;
	next if (/^>/);
	my @seq = split( '', $_);
	for ( @seq ) {
		$data{$_}++;
	}
}

$data{'total'} = sum( map { $data{$_} } keys %data );
$data{'GC'} = $data{'G'} + $data{'C'};

say $data{'GC'}/$data{'total'};

计算所得文件之后结果为:0.51911738950782

3、谷歌一下可知人类染色体GC含量为
#Sequence   GC content
chr1          0.43
chr2          0.40
chr3          0.40
chr4          0.38
chr5          0.40
chr6          0.40
chr7          0.41
chr8          0.40
chr9          0.43
chr10         0.42
chr11         0.42
chr12         0.41
chr13         0.40
chr14         0.43
chr15         0.44
chr16         0.45
chr17         0.46
chr18         0.40
chr19         0.48
chr20         0.44
chr21         0.43
chr22         0.49
chrX          0.40
chrY          0.46
chrM          0.44

出处https://www.biostars.org/p/16169/

可知平均GC含量在45%左右,可见CDS的CG含量确实高于平均水平。

这就是我分析的全过程,欢迎指教。
回复 支持 4 反对 0

使用道具 举报

0

主题

15

帖子

72

积分

注册会员

Rank: 2

积分
72
发表于 2017-2-11 00:49:14 | 显示全部楼层
考虑的不够细,大致统计了下。。

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

use strict;
use warnings;

my ($genome,$gtf);

if(@ARGV < 2){
  die "Usage: perl $0 <Genome> <GTF>";	
}else{
  ($genome,$gtf) = @ARGV;	
}

open GM, '<'.$genome;
open OUT, '>'."$genome.GCstat";


my $tit;
my %seq;

while(<GM>){
  chomp;
  if(/>(\S+)/){
    $tit = $1;
  }else{
    $seq{$tit} .= $_;	
  }
}
close GM;

open GTF, $gtf;
#my ($chrom,$element,$start,$end);
my %exon_seq;
while(<GTF>){
  chomp;
  next if /^\#/;
  my ($chrom,$element,$start,$end) = (split/\t/)[0,2,3,4];
  if($element eq 'exon'){
    $exon_seq{$chrom} .= substr($seq{$chrom},($start - 1),($end - $start + 1));
  }
}

print OUT "Chrom\tAll_GC\tAll_GC_per\tExon_GC\tExon_GC_per\n";
foreach my $chrom (sort keys%seq){
  my $chr_seq = $seq{$chrom};
  my $all_exon_seq = $exon_seq{$chrom};
  print OUT $chrom."\t".stat_GC($chr_seq)."\t".stat_GC($all_exon_seq)."\n";
}


sub stat_GC {
  my ($in_seq) = @_;
  my %basecount;
  my ($GC_count,$GC_percent);
  my $seq_len = length $in_seq;
  map {$basecount{uc($_)}++} (split//,$in_seq);
  $GC_count = $basecount{'G'} + $basecount{'C'};
  if(defined $basecount{'N'}){
    $GC_percent = $GC_count/($seq_len - $basecount{'N'});
  }else{
    $GC_percent = $GC_count/$seq_len;	
  }
  return($GC_count."\t".$GC_percent);
}


本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-3-27 12:52 , Processed in 0.026681 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.