最新更新：
最簡(jiǎn)單的方法還是R 一句話d <-data1[,data2$Header] （data1就是這里的1.txt, data2就是2.txt）
今天遇到了一個(gè)小問(wèn)題颖变，想把轉(zhuǎn)錄組read count矩陣中指定樣品（指定列）的表達(dá)量挑選出來(lái)听想，總共從500多個(gè)樣中選200個(gè)马胧，數(shù)據(jù)量在500Mb左右佩脊。為了簡(jiǎn)化問(wèn)題，我先測(cè)試了一下威彰。
假設(shè)有1.txt和2.txt兩個(gè)文件歇盼，格式如下：
$ 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
現(xiàn)在根據(jù)2.txt里指定列的信息從1.txt里挑第一列和第二列出來(lái)，最終想得到這樣的結(jié)果：
1 2
a b
g h
方案1. shell腳本（冗雜且提取失敳础）
我寫(xiě)了一個(gè)簡(jiǎn)易的shell腳本可惜不成功
$ more test.sh
a=`awk '{print NF}' 1.txt` #統(tǒng)計(jì)1.txt的列數(shù)
b=`wc -l |2.txt` #統(tǒng)計(jì)2.txt的行數(shù)
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'}'` #逐個(gè)讀取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結(jié)果不太對(duì)，提取列之后要paste 還是不方便組合
琢磨了一陣后侍匙，我向生信技能樹(shù)的小伙伴們求助叮雳，果然群體的智慧是無(wú)窮的~
方案2. awk提攘辈弧（感謝王詩(shī)翔的建議和幫助）
原話：
linux處理文本的核心是以行為基礎(chǔ)，我的意思是利用現(xiàn)有的腳本將列變成行唬滑，然后使用join拼接
i:
i:$i}END{for(i=0;i++<NF;)print a[i]}' 行列轉(zhuǎn)換的命令網(wǎng)上就有棺弊，也可以自己寫(xiě)
$ 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
不得不說(shuō)學(xué)好這里shell命令真的方便，join實(shí)現(xiàn)的功能之前我是用python腳本弄的......
不過(guò)尊勿，這里有個(gè)問(wèn)題需要注意畜侦，join是按行提取的旋膳，如果有一行在1.txt和2.txt里面不匹配，就會(huì)停止檢索验懊。
比如义图，2.txt里面多加一行9（1.txt里面沒(méi)有）
$ more 2.txt
1
2
9
3
$join 1.txt 2.txt
1 a g
2 b h
后面的第3行3 c i就沒(méi)有被提取出來(lái)
方案3. R包dplyr select()提取（感謝嚴(yán)濤的建議和幫助娃承，這是他的個(gè)人R學(xué)習(xí)筆記里的一部分）
首先在2.txt首行加個(gè)Header 方便提取
$ more 2.txt
Header
1
2
$ R
>library(dplyr)
>3.txt <- 1.txt %>% select(one_of(dput(as.character(2.txt$Header))))
這里%>% 是管道函數(shù)，把左邊文件的值發(fā)送給右邊文件桶蛔，并作為右邊文件件表達(dá)式的第一個(gè)參數(shù), select()允許我們快速通過(guò)變量名對(duì)數(shù)據(jù)集取子集漫谷，后面的看的不是很懂
推薦一篇王詩(shī)翔寫(xiě)的介紹dplyr的博客詳細(xì)了解一下
使用dplyr進(jìn)行數(shù)據(jù)轉(zhuǎn)換
方案4.python按行提忍蚴尽（我之前用的腳本）
還是要先轉(zhuǎn)置，然后再提,不過(guò)只提了前兩列
#!/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腳本處理的思路有很多俺祠，這里是用先轉(zhuǎn)置成行再存數(shù)組匹配，大致是這樣
轉(zhuǎn)置
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;
}
方案6. plink提饶韪住（樣品很多的終極選擇）
plink --vcf input.vcf --noweb --keep-fam id.txt --recode --make-be
d --out output
id.txt里面按行寫(xiě)入要提取的列名就行，很方便
網(wǎng)盤(pán)下載鏈接 鏈接:
Extract.pl
https://pan.baidu.com/s/1pCdM0tch8Jl2QeZO4clQ9g 密碼:97hy
Trans.pl
https://pan.baidu.com/s/18aMepxBrk7EHbAXmTfVyOw 密碼:b10o
非常感謝給予幫助的小伙伴吐葱，當(dāng)然還有很多方法可以實(shí)現(xiàn)校翔。P.S.最初的簡(jiǎn)陋版shell如果有高手看出問(wèn)題防症，歡迎留言給我~~
