用于查找文件中最长和最短的单词或字符串的脚本？

Question

我正在做基因组学。我有一个包含 FASTA 格式读取的文件。这些是基因。每个基因称为读取或重叠群。每个重叠群以标题开头，后跟特定长度的字母或核苷酸，例如：ACTG。我想确定该文件中最长的重叠群和最短的重叠群或读取或基因。请告诉我一个 ubuntu 脚本来查找此类重叠群。每个重叠群或读取都采用这种 FASTA 格式，如下所示：

>Locus_1000_Transcript_1/1_Confidence_0.000_Length_648 FTBs=645 (Header)
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCAgacGcccGttCacAAaAaC

标题具有定义的格式，显示基因位点和基因数量，并且每个重叠群或读段之间会有空格。文件中的每个读段或重叠群都将以与上述类型相同的标题开头，但值可能不同。每个重叠群或读段都以 > 符号开头。可能会有长度相同的重叠群。– Science 3 分钟前

Answer 1

假设LengthFASTA 标题中的值是正确的，我会从那里提取它们：

sed -nre 's/^>.*_Length_([0-9]+) .*/\1/p' \

然后按数字排序

| sort -n \

然后输出第一行和最后一行

| sed -ne '1p;$p'

有一条声明：

sed -nre 's/^>.*Length_([0-9]+) .*/\1/p' | sort -n | sed -ne '1p;$p'

---

如果标题中声明的长度不可信，那么为了计算 FASTA 序列的长度，我首先将它们转换为联合国，然后将每隔一行的行长度打印到与sort | sed上面相同的过滤器中：

uf | awk 'NR%2==0 {print length}' | sort -n | sed -n '1p;$p'

在哪里uf可以找到简单的 bash 脚本这里。

---

注意：这两个单行命令都是过滤器，也就是说它们从标准输入读取输入并写入标准输出。使用cat它们来提供文件（或wget -O -从互联网上提供文件）。

Answer 2

这个 python 脚本建立了记录字典，并使用线性搜索来查找文件中哪个最短、哪个最长。请注意，这会忽略两个具有相同值的 contig 的情况（尽管也可以实现）。

代码：

#!/usr/bin/env python3

import sys
def main():
    records = {}
    current_length = 0
    current_contig = ''
    with open(sys.argv[1]) as f:
        for index,line in enumerate(f,1):
            if line == '\n': continue
            if line.startswith('>'):
               if current_contig != line:
                   records[current_contig] = current_length
                   current_contig = line.strip()
                   current_length = 0
            else:
               current_length = current_length + len(line.strip())
    records[current_contig] = current_length
    records.pop('')

    shortest_contig = None
    longest_contig = None 
    longest_val = 0
    shortest_val = float("inf")
    for contig,length in records.items(): 
        if length < shortest_val:
            shortest_val = length
            shortest_contig = contig
        if length > longest_val:
            longest_val = length
            longest_contig = contig
    print('Longest: ' + longest_contig)
    print('Shortest: ' + shortest_contig)


if __name__ == '__main__': main()

测试运行：

$ cat input.txt                                                                                                                  
> Entry 1
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCAgacGcccGttCacAAaAaC

> Entry 2
CcGccttggtaacctCgccAGcatATtgagcTttGGatccGGaTggtcgtaGaAtgGCaaG
GcaGgagAgAgtgtctaatgtggCgccGctctgtAccCgGgGGgTAACaAtgAATTtGCga
CgaCGtggTAtGcCcttCGttgAaacccTtaTtagttGgAGCcGctAtgtggcgGTccaat
TaTcaagtAttTcCCACaTcttgAagCgcttcTgGATgTacgCatactatgggTtgacgtt
AGtGtAgCcgAgattTCacaGtAgctcCGAACGgtgGTagCA
$ python3 contigs.py  input.txt                                                                                                   
Longest: > Entry 1
Shortest: > Entry 2