我遇到了尝试在unix中加入两个数据集的问题,可以使用你的帮助。我花了很长时间在论坛上寻找解决方案,但空手而归。
我在一个数据集中有一个入藏号列表,需要将它们转换为基因符号。为了做到这一点,我下载了gene2accession.gz from NCBI。未压缩的文件是〜7Gb,所以首先我从这个数据集中删除了加入和基因符号
cut -f 2,16 gene2accession > accession2genesymbol
根据wc -l accession2genesymbol
有大约7000万行,有许多重复,所以我用sort accession2genesymbol | uniq
删除了这些,导致了大约2000万行。
现在通常我会在R中做一个inner_join()
using the dplyr
package(从x返回y中匹配值的所有行,以及x和y中的所有列);但是,这个数据集对我来说太大了。
以下是未分类的accession2genesymbol数据集的示例:
100000492 mafaa
1000004 XCC3444
110047139 LOC110047139
110047140 LOC110047140
9951915 LOAG_14435
9951916 LOAG_14436
999999 gndA
999 CDH1
9 NAT1
未分类的Accessions的简短示例如下所示(对于整个数据集-1,576行see the gist):
Accessions
10047140
100913206
10092617
10190704
10190704
103471987
103471997
103472005
103472005
105990514
45006951
45006986
45006986
45007007
45007007
4501883
4501887
94721250
94721261
9558733
9845516
98986457
98986457
98986464
99028871
9910242
9951915
9966805
9966827
9966867
9994185
编辑:只有加入9951915和110047140这里有匹配所以我的预期输出将是:
9951915 LOAG_14435
110047140 LOC110047140
由于没有使用unix进行数据操作/加入,我搜索了Stack Overflow的类似问题。
例如this one。我的理解是unix join
函数只能在文件排序时使用,所以我尝试了以下方法:
join -t "\t" <(dos2unix <accession) <(dos2unix <accession2genesymbol.txt)
也许这不起作用,因为我在两个数据集中都需要完全相同的行号(即如果数据集的第100行与数据集2的行100不匹配则不行)但也许我错了,这对某些人不起作用其他原因?
也许awk
是一个更好的解决方案,所以我尝试了this post的建议:
awk '{a[$1]=a[$1] FS $2} END {for (i in a) print i a[i]}' accession accession2genesymbol | sort > file3
这会生成一个包含大约2000万行的文件,因为我的加入只有9000行,所以我预计会有9000行(例如,如果这些访问不再存在,则可能更少)。
我尝试了第一篇文章中的另一个awk
解决方案:
awk -F, 'FNR==NR{a[$1];next}($1 in a){print $2}' accession accession2genesymbol > file3
awk: warning: escape sequence `\+' treated as plain `+'
但我最终得到一个空文件。
我很欣赏awk(病房)解决方案,python或任何可以帮助我解决这个问题的方法。非常感谢你。
我们还没有看到你的原始gene2accession
文件的样本,但让我们假设它是一个制表符分隔的字段,第二列中有accession
,第16列中的gene
(因为这是你的cut
正在选择的)带有标题行。让我们假设您的Accessions
文件并非绝对巨大。
鉴于此,这个脚本应该做你想要的:
awk -F'\t' 'NR==FNR{a[$1];next} ($2 in a) && !seen[$2]++{print $2, $16}' Accessions gene2accession
但你可以试试看它是否更快:
awk -F'\t' 'NR==FNR{a[$1];next} $2 in a{print $2, $16}' Accessions <(sort -u -t'\t' -k2,2 gene2accession)
如果是,你想要一个sort
输出的中间文件用于后续运行:
sort -u -t'\t' -k2,2 gene2accession > unq_gene2accession &&
awk -F'\t' 'NR==FNR{a[$1];next} $2 in a{print $2, $16}' Accessions unq_gene2accession
join
应该适合你的情况。由于您的输入文件没有匹配,因此这是一个组成示例并使用您的地图文件
$ head file
100000009
100000061
100000030
$ join <(sed 1d map) <(sort file)
100000009 sema5bb+
100000030 btr24+
100000061 si:ch211-133n4.9+
假设您的map
文件已经排序,您需要删除标题sed 1d
并需要对您的输入file
进行排序。请注意,排序应该是数字或词汇。
另一种不需要分类的替代方案是使用grep
$ grep -wFf file map
100000009 sema5bb+
100000030 btr24+
100000061 si:ch211-133n4.9+
如果数字和代码的格式不同,则不会出现错误匹配。