一、对每一个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