折腾了几周,基于perl编程,根据基因长度及其在染色体上的位置分布信息,在全基因组上我识别出了42641个组蛋白修饰模式组合(也就是某些人说的组蛋白密码),然后去除冗余,得到18506个不重复的组蛋白修饰模式。统计分析发现,其中大部分模式在全基因组上只出现了1次。出现次数大于等于4次的模式一共有927个。之前我在全基因组上用5000bp的半重叠窗去采集24条染色体上的组蛋白修饰模式,共得到854个模式。两者相比很多模式都比较接近。说明这些模式的确在染色体上发生了。下一步我想研究一下这些组蛋白修饰模式与DNA甲基化之间的关系。DNA甲基化以前没怎么接触过,希望好运~~,不说了,贴出批处理24条染色体,39种组蛋白修饰的perl代码。
#!/bin/perl
use strict;
my $file_name;
my @doc_name=qw /chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 chr20 chr21 chr22 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrM chrX chrY /;
foreach $file_name (@doc_name) {
my %Data;
open(INFILE,"$file_name.txt");
while(<INFILE>) {
chomp;
if (/^\w+\t(\d+)\t(\w+)\t(\d+)\t(\d+)/) {
my $gene = "$2-$1($3-$4)"; #Chromsome-geneID(Start-End)
push @{$Data{$gene}},$3,$4;
}
}
close(INFILE);
open(DATA,"$file_name");
my @arr;
while(<DATA>)
{
push @arr,[split];
}
close DATA;
my @histone=qw /CD4-H2BK120ac CD4-H4K91ac CD4-H2BK20ac CD4-H3K4ac CD4-H3K36ac CD4-H2BK5ac CD4-H3K27ac CD4-H3K18ac CD4-H3K9ac CD4-H2BK12ac H3K4me3 H2AZ CD4-H4K5ac CD4-H4K8ac H3K9me1 H3K4me2 H3K4me1 H2BK5me1 H4K20me1 H3K79me3 H3K79me2 H3K79me1 CD4-H4K12ac CD4-H4K16ac CD4-H2AK5ac H4K20me3 H4R3me2 H3R2me2 H3R2me1 H3K27me1 H3K36me3 CD4-H3K23ac CD4-H2AK9ac CD4-H3K14ac H3K27me3 H3K36me1 H3K27me2 H3K9me3 H3K9me2/;
my %hash;
my %index;
my @temp;
for (0..38) {
$hash{$histone[$_]}=1;
$index{$histone[$_]}=$_;
$temp[$_]=0;
}
open(OUT,">out_$file_name");
foreach my $Locu (keys %Data) {
my $spaceflag=0; #控制多余空行的输出
for(my $i=0;$i<=$#arr;$i++){
if($arr[$i]->[1] >= $Data{$Locu}->[0] && $arr[$i]->[1]<=$Data{$Locu}->[1] ){
#print OUT $Locu."\t"; #输出基因名及其位点信息
$temp[$index{$arr[$i]->[0]}]=$hash{$arr[$i]->[0]}; #$index{$arr[$i]->[0]}取值0~38;$hash{$arr[$i]->[0]}取值39个修饰
$spaceflag=1;
}
}
if($spaceflag==1){
print OUT $Locu."t"; #统计重复时此行不要。
print OUT join("t",@temp);
print OUT "n";
}
for (0..38) { #每循环一次,变量@temp归零。
$temp[$_]=0;
}
}
close OUT;
}
https://blog.sciencenet.cn/blog-634847-558844.html
上一篇:
perldoc真是一个好东西~