今天遇到了一个小问题，想把转录组read count矩阵中指定样品（指定列）的表达量挑选出来，总共从500多个样中选200个，数据量在500Mb左右。为了简化问题，我先测试了一下。
假设有1.txt和2.txt两个文件，格式如下：
$ more 1.txt
1 2 3 4 5 6
a b c d e f
g h i j k l
$ more 2.txt
1
2
现在根据2.txt里指定列的信息从1.txt里挑第一列和第二列出来，最终想得到这样的结果：
1 2
a b
g h
方案1. shell脚本（冗杂且提取失败）
我写了一个简易的shell脚本可惜不成功
$ more test.sh
a=`awk '{print NF}' 1.txt` #统计1.txt的列数
b=`wc -l |2.txt` #统计2.txt的行数
for (( i=1;i<=$a;i++))
do
for(( j=1;j<=$b;j++ ))
do
h=`cat 1.txt |awk 'NR==1{print}'|awk '{print '$i'}'` #逐个读取1.txt第一列
k=`cat 2.txt |awk 'NR=='$j'{print}'` #读取2.txt的每一行
if [[ h -eq k ]];
# then echo $k
then echo `cat 1.txt |awk '{print '$i'}'`
fi
done
$ sh test.sh
1 1
2 2
echo结果不太对，提取列之后要paste 还是不方便组合
琢磨了一阵后，我向生信技能树的小伙伴们求助，果然群体的智慧是无穷的~
方案2. awk提取（感谢王诗翔的建议和帮助）
原话：
linux处理文本的核心是以行为基础，我的意思是利用现有的脚本将列变成行，然后使用join拼接 $ awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}' 1.txt | join - 2.txt | awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}' 行列转换的命令网上就有，也可以自己写
$ awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}' 1.txt | join - 2.txt | awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}'
1 2
a b
g h
拆解一下 先把行和列置换，然后再用join命令按行匹配,再置换一次就好了
$ awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}' 1.txt >01
1 a g
2 b h
3 c i
4 d j
5 e k
6 f l
$ join 01 2.txt >02
1 a g
2 b h
$ awk '{for(i=0;++i<=NF;)a[i]=a[i]?a[i] FS $i:$i}END{for(i=0;i++<NF;)print a[i]}' 02
1 2
a b
g h
不得不说学好这里shell命令真的方便，join实现的功能之前我是用python脚本弄的......
不过，这里有个问题需要注意，join是按行提取的，如果有一行在1.txt和2.txt里面不匹配，就会停止检索。
比如，2.txt里面多加一行9（1.txt里面没有）
$ more 2.txt
1
2
9
3
$join 1.txt 2.txt
1 a g
2 b h
后面的第3行3 c i就没有被提取出来
方案3. R包dplyr select()提取（感谢严涛的建议和帮助，这是他的个人R学习笔记里的一部分）
首先在2.txt首行加个Header 方便提取
$ more 2.txt
Header
1
2
$ R
>library(dplyr)
>3.txt <- 1.txt %>% select(one_of(dput(as.character(2.txt$Header))))
这里%>% 是管道函数，把左边文件的值发送给右边文件，并作为右边文件件表达式的第一个参数, select()允许我们快速通过变量名对数据集取子集，后面的看的不是很懂
推荐一篇王诗翔写的介绍dplyr的博客详细了解一下
使用dplyr进行数据转换
方案4.python按行提取（我之前用的脚本）
还是要先转置，然后再提,不过只提了前两列
#!/usr/bin/python
file1=open("",'r')
file2=open("1.txt",'r')
file3=open"2.txt",'w')
file1_dict={}
while 1:
line1=file1.readline()
if not line1:
break
lin=line1.strip('\n')
lin1=lin.split('\t')
file1_dict[lin1[0]]=lin1[1]
while 1:
line2=file2.readline()
if not line2:
break
line=line2.strip('\n')
if line in file1_dict:
value=file1_dict[line]
file3.write(line+'\t'+value+'\n')
file1.close()
file2.close()
file3.close()
方案5. perl脚本提取（感谢刘帅的建议和帮助）
perl脚本处理的思路有很多，这里是用先转置成行再存数组匹配，大致是这样
转置
while (my $tem=<IN>){
chomp $tem;
my @ll=split /\t/,$tem;
push @sample,$ll[0];
for my $i (1..$#ll){
push @{$snp{$snp_name[$i]}},$ll[$i];
}
}
提取行
while (<IN1>) {
chomp;
my @a=split/\s+/,$_;
push @sample,$a[0];
}
while (<IN2>) {
chomp;
my @b=split/\s+/,$_;
foreach $i(@sample) {
if ($i eq $b[0]) {
print OUT "$_\n";
}
完整版看这里：
Trans.pl
my @id;
my @chr;
my @pos;
my $head1=<IN>;
chomp $head1;
my @snp_name=split /\t/,$head1;
while (my $tem=<IN>){
chomp $tem;
my @ll=split /\t/,$tem;
push @sample,$ll[0];
for my $i (0..$#ll){
push @{$snp{$snp_name[$i]}},$ll[$i];
}
}
open (OUT,">$outfile") || die "Can't creat $outfile, $!\n" ;;
for my $i (0..$#snp_name) {
my $content=join("\t",@{$snp{$snp_name[$i]}});
print OUT "$snp_name[$i]\t",$content,"\n";
}
sub USAGE {#
my $usage=<<"USAGE";
ProgramName: Transpose of Matrix
Version: $version
Contact: Shuai Liu <ls2106\@msstate.edu>;
Program Date: 2018.6.9
Usage:
Options:
-infile <file> input file,forced
-outfile <file> output file,forced
-h Help
USAGE
print $usage;
exit;
}
Extract.pl
#!/usr/bin/perl -w
use strict;
use warnings;
use Getopt::Long;
my $version="1.0";
#######################################################################################
my ($list,$infile,$outfile);
GetOptions(
"help|?" =>\&USAGE,
"list:s"=>\$list,
"infile:s"=>\$infile,
"outfile:s"=>\$outfile,
) or &USAGE;
&USAGE unless ($list||$infile||$outfile);
######################### vcffilter & vcf imputation ###############################
open (IN1, "<$list") || die "Can't creat $list, $!\n" ;
my @sample;
my $a;
my $i;
my @b;
my @a;
while (<IN1>) {
chomp;
my @a=split/\s+/,$_;
push @sample,$a[0];
}
open (IN2, "<$infile") || die "Can't creat $list, $!\n" ;
open (OUT, ">$outfile") || die "Can't creat $list, $!\n" ;
while (<IN2>) {
chomp;
my @b=split/\s+/,$_;
foreach $i(@sample) {
if ($i eq $b[0]) {
print OUT "$_\n";
}
}
}
sub USAGE {#
my $usage=<<"USAGE";
ProgramName: Extract rows from list
Version: $version
Contact: Shuai Liu <ls2106\@msstate.edu>;
Program Date: 2018.6.9
Usage:
Options:
-list <file> list file,forced
-infile <file> input file,forced
-outfile <file> output file,forced
-h Help
USAGE
print $usage;
exit;
}
网盘下载链接 链接:
Extract.pl
https://pan.baidu.com/s/1pCdM0tch8Jl2QeZO4clQ9g 密码:97hy
Trans.pl
https://pan.baidu.com/s/18aMepxBrk7EHbAXmTfVyOw 密码:b10o
非常感谢给予帮助的小伙伴，当然还有很多方法可以实现。P.S.最初的简陋版shell如果有高手看出问题，欢迎留言给我~~
网友评论