标签：22 23 my perl subData flag line SVA open

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

标签：22,23,my,perl,subData,flag,line,SVA,open
来源： https://www.cnblogs.com/zjuer/p/16513644.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。