一、对每一个repeat family的转录因子的富集进行提取。
由于我好久没有用perl了，每次想要对文本文件进行提取的时候，总会想到perl。但是由于对代码的生疏（好久没有用了），主要的意思是明白的，但是具体的指令写起来就特别的费劲。
现在开始一步一步的写，给自己一个小时的时间，写到晚上十点。
实现的目标是：（1）输出提取后的ID名；
（2）每一个序列区间，下面有一堆的ID名。
以下使用极其蹩脚的perl代码解决了。
open (MARK, "< human_specific_repeat_uniq_intervals.result") or die "can not open it!";
open (COUNT, "> count_result.txt") or die "can not open it!";
$/="\n";
my $flag;
#my %data;
while (<MARK>){
$line =$_;
$flag = substr($line,0,2);
if ($flag eq "##"){ #bug!在这里
#@repeat = split /\s+/,$line;
#$data{$repeat[5]} = "";
chomp $line;
print COUNT $line,"\n";
}
else
{
chomp $line;
@list = split(/\s+/,$line);
my @T = split /_/,$list[10];
my @TF = split /\//,$T[2];
print COUNT $TF[1],"\n";
}
}
close(MARK);
close(COUNT);
怎么可能一下子找到正确答案的你说，我们总要经历这样一个试错的过程。不断地化繁去简，接近真相。
##chr10 100092148 100093575 SVA_F SVA SVA
32850
32853
32856
37635
41187
43676
43684
48334
48930
49606
64738
64739
67858
70095
73982
二、转录因子的值的转换。
接下来，对上述文件进行读取，加上转录因子的那个文件，对其进行替换。
> data<-read.table("human_factor_full_QC.txt",sep="\t",fill=T,header=T,na.strings="",comment.char="#",quote="")
> subData<-subset(data,select=c("DCid","Factor"))
> head(subData
+ )
DCid Factor
1 1 BTAF1
2 2 GAPDH
3 4 EGR1
4 6 TCF4
5 8 TCF4
6 9 TCF4
> write.table(subData,"DCid_TF.txt",quote=F,row.names=F)
然后想办法把数字替换成转录因子的名字，是不是我们计数后再提取会更好操作一些呢？
这里的处理则相对比较简单，使用merge函数即可满足要求。
> data<-read.table("SVA_result_v1.txt")
> head(data)
V1
1 1
2 1006
3 1007
4 1010
5 1012
6 1082
> TF<-read.table("DCid_TF.txt",sep="\t",header=T)
> head(TF)
DCid.Factor
1 1 BTAF1
2 2 GAPDH
3 4 EGR1
4 6 TCF4
5 8 TCF4
6 9 TCF4
> colnames(data)<-"DCid"
> subData<-merge(data,TF)
> head(subData)
DCid Factor
1 1 BTAF1
2 2 GAPDH
3 4 EGR1
4 8 TCF4
5 9 TCF4
6 11 TCF4
> dim(subData)
[1] 4908 2
> write.csv(subData,"SVA_TF.csv",row.names=F
其他几个文件同样的处理。
到这里我们的需求基本上就结束了，问题也解决了。我们接下来就想说，对这些转录因子进行功能注释。（明天的工作内容）
三、汇总同一家族的转录因子的计数的情况。
之前使用的是python的字典来存储这种数据格式，我总觉得perl中的哈希也是同样的数据格式。这次使用哈希来试一试。
蹩脚小张终于把代码写出来了，真的是边玩边写的。这里主要是两个思路解决了我的问题。（1）首先，想到的是以列表的形式先把值拿出来（尽管存在重复）；（2）其次，使用了Data::Dumper这个库，对哈希进行规范化的输出；
#!perl
use Data::Dumper; #我对一些库掌握的还很少；现在只会一些基础的语法；
open (INPUT, "< count_result.txt") or die "can not open it!";
open (COUNT, "> SVA_result.txt") or die "can not open it!";
$/="##";
my @flag;
#创建一个已知的字典:LINE/SINE/LTR/SVA
my %count=("LINE"=>[],"SINE"=>[],"LTR"=>[],"SVA"=>[]);
readline <INPUT>;
while (<INPUT>){
$line =$_;
chomp $line;
@flag = split /\n+/,$line;; #提取该行数据
@repeat = split /\s+/,$flag[0];
my $len=@flag;
for( my $a = 1; $a < $len; $a = $a + 1 ){
push @{$count{$repeat[5]}}, $flag[$a];
}
}
# while(my ($key, $val) = each(%count)) {
# print $key,$val[1];
# }
# @output = $count{"SVA"};
# print @output;
print COUNT Dumper(\$count{"SVA"});
现在对结果数据进行整理和计数。
首先使用文本编辑器，对数据进行初步的处理，然后去除重复。
cat LTR_result.txt |sort |uniq >LTR_result_v1.txt
