今天在做miRNA分析，被一个问题困住了
其中有一个步骤是用CD-hit把序列聚类去冗余，得到聚类的文件和去冗余后的文件。其中聚类的文件是这样的：
>Cluster 0
0 35aa, >NS500576:848:HJ5LLBGXH:1:11101:2696:1553... *
>Cluster 1
0 34aa, >NS500576:848:HJ5LLBGXH:1:11101:19894:1287... at 100.00%
1 35aa, >NS500576:848:HJ5LLBGXH:1:11101:5101:1785... *
2 34aa, >NS500576:848:HJ5LLBGXH:1:11101:15903:2331... at 100.00%
3 34aa, >NS500576:848:HJ5LLBGXH:1:11101:24910:2474... at 100.00%
4 35aa, >NS500576:848:HJ5LLBGXH:1:11101:7539:2794... at 100.00%
5 35aa, >NS500576:848:HJ5LLBGXH:1:11101:18761:3213... at 100.00%
6 35aa, >NS500576:848:HJ5LLBGXH:1:11101:13988:3650... at 100.00%
7 34aa, >NS500576:848:HJ5LLBGXH:1:11101:24562:3668... at 100.00%
8 34aa, >NS500576:848:HJ5LLBGXH:1:11101:8305:3802... at 100.00%
9 35aa, >NS500576:848:HJ5LLBGXH:1:11101:22674:3841... at 100.00%
而我想要的是什么结果呢？
我想要的是每一个类的最长序列对应这个类的序列数目，也就是类似下面的内容:
NS500576:848:HJ5LLBGXH:1:11101:5101:1785 100
NS500576:848:HJ5LLBGXH:1:11101:2696:1553 1
而我最熟悉的就是shell命令，怎么用shell不借助其他语言就能生成我想要的文件呢？
在此之前查看了非冗余的序列文件，发现它的序列的顺序并不是和聚类的是一致的。
左思右想.......................
突然：灵感来了！
每个聚类的最长序列都有一个特点，他们最后都是以*结尾的，那这样我就可以把最长序列和每个聚类对应起来了。
然后怎么获取一个聚类有多少条序列呢？（这是第二个问题）
哎！仔细看这个文件，每一个聚类名字前面一行对应的是上一个聚类的最后一个序列的名字以及编号，其编号加上1就是上一个聚类有多少条序列。所以，那我只需要提取每个聚类名的上一行（注意第一行的上一行是空的），把其对应的编号加上1，再和我第一步提取出来的带*号的序列合并，去掉没用的字符，好像是可以生成我想要的文件哎。
实践时间！
- 第一步:提取每个聚类对应的最长序列
grep "*" V0-1.100.clstr|wc -l #先查看一下有多少个序列
361700
wc -l V0-1.100.fa #查看CD-hit输出的非冗余序列的行数
723400 #刚好是上面的两倍，因为该序列包括序列名称以及序列的具体信息，所以证明我们按*来提取序列名称这个方法没有问题
grep "*" V0-1.100.clstr|tr ">" "\t"|cut -f3|sed "s/\.\.\.\ \*//g">V0-1.100.cls.longseq #把最长序列提取出来并去掉没用的字符
- 第二步：提取每个聚类上一行（主要是为了获取编号）
sed -n '/Cluster/{x;p};h' V0-1.100.clstr |cut -f1|sed 1d|awk '{print$1+1}'>V0-1.100.cls.number
wc -l V0-1.100.cls.number
361699 #发现少一个，应该是最后一个聚类没有取到
tail -n1 V0-1.100.clstr
0 17aa, >NS500576:848:HJ5LLBGXH:4:23612:24745:20134... *
echo "1" >>V0-1.100.cls.number #把最后一个聚类的数目加到前面得到的文件上面
wc -l V0-1.100.cls.number
361700 #再检查一遍 结果是对的
- 第三步：合并前两步得到的文件
paste V0-1.100.cls.longseq V0-1.100.cls.number >V0-1.100.cls.longseq.number
- 最后得到了我想要的文件
less -S V0-1.100.cls.longseq.number
NS500576:848:HJ5LLBGXH:1:11101:2696:1553 1
NS500576:848:HJ5LLBGXH:1:11101:5101:1785 22885
NS500576:848:HJ5LLBGXH:1:11101:6861:1893 293
NS500576:848:HJ5LLBGXH:1:11101:22092:2106 2
NS500576:848:HJ5LLBGXH:1:11101:12394:2309 433
NS500576:848:HJ5LLBGXH:1:11101:18059:2391 1487
NS500576:848:HJ5LLBGXH:1:11101:13072:3057 2216
NS500576:848:HJ5LLBGXH:1:11101:14749:3111 76957
NS500576:848:HJ5LLBGXH:1:11101:7715:3538 462
NS500576:848:HJ5LLBGXH:1:11101:5549:3540 182
NS500576:848:HJ5LLBGXH:1:11101:19431:3570 115
NS500576:848:HJ5LLBGXH:1:11101:20139:3678 775
NS500576:848:HJ5LLBGXH:1:11101:21482:3830 1142
NS500576:848:HJ5LLBGXH:1:11101:7965:4295 678
NS500576:848:HJ5LLBGXH:1:11101:25840:4929 1
想说一句，我的智商还可以！
网友评论