主页 / unix / 问题

问题[bioinformatics](unix)

Martin Hope
compbiostats
Asked: 2022-09-22 05:19:40 +0800 CST
  • 2

我是新手GNU sed（在 macOS 上运行）并希望缩短一些文件头（~50 K 头，78.3 Mb）。

我正在尝试通过指定带有-i标志的备份扩展名来修改 FASTA 文件。

到目前为止我有

sed -i.bak 's/^([^|]+).[^|]+(.*)/\1\2/' file.fas

这应该创建一个名为file.fas.bak.

但是，我收到了错误

sed: 1: "s/^([^|]+).[^|]+(.*)/\1\2/": \1 not defined in the RE

注意sed -re 's/^([^|]+).[^|]+(.*)/\1\2/' file.fas正确打印到屏幕。

在这种情况下如何打印到备份文件有什么想法吗？

sed bioinformatics
Martin Hope
geneteics_diva
Asked: 2022-07-07 22:14:43 +0800 CST
  • 1

我是脚本新手，所以我需要帮助。

我正在运行一个基于四个独立测试（--max-maf）的程序，该测试需要输入文件名和输出文件名。下面只是对该程序的一般描述。文件名是我在下一个代码块中详细描述的输入。

epacts group --groupf filename.grp --vcf filename.vcf.gz --ped filename.ped --max-maf 0.05 --kin filename --test emmaxCMC --out BcA/filename-CMC-0.05
epacts group --groupf filename.grp --vcf filename.vcf.gz --ped filename.ped --max-maf 0.03 --kin filename --test emmaxCMC --out BcA/filename-CMC-0.03
epacts group --groupf filename.grp --vcf filename.vcf.gz --ped filename.ped --max-maf 0.02 --kin filename --test emmaxCMC --out BcA/filename-CMC-0.02
epacts group --groupf filename.grp --vcf filename.vcf.gz --ped filename.ped --max-maf 0.01 --kin filename --test emmaxCMC --out BcA/filename-CMC-0.01

我在同一个目录中有几个具有相同基本名称的输入文件，但一个文件 (CDES_MyopV1.ped) 具有相同的基本名称，但后面有一个唯一标识符。该文件将在 --ped 命令之后执行 (--ped CDES_MyopV1.ped)

CDES-genes.grp 
CDES.vcf.gz 
CDES_MyopV1.ped 
CDES.kinf

我尝试通过根据基本名称“CDES”查找上面列出的输入文件来创建一个执行程序的脚本，但是我意识到我需要 .ped 文件来包含基本名称和其后的唯一标识符（CDES_MyopV1）此外，对于每个输出文件，我希望将 CDES_MyopV1 连接到输出文件名。

这是我迄今为止尝试过的：

declare -a files=("CDES")

for element in ${files[@]}
do
   epacts group --groupf $element-genes.grp --vcf $element.vcf.gz --ped $element.ped --max-maf 0.05 --kin $element.kinf  --test emmaxCMC --out BcA/$element-CMC-0.05
   epacts group --groupf $element-genes.grp --vcf $element.vcf.gz --ped $element.ped --max-maf 0.03 --kin $element.kinf  --test emmaxCMC --out BcA/$element-CMC-0.01 
   epacts group --groupf $element-genes.grp --vcf $element.vcf.gz --ped $element.ped --max-maf 0.02 --kin $element.kinf --test emmaxCMC --out BcA/$element-CMC-0.05 
   epacts group --groupf $element-genes.grp --vcf $element.vcf.gz --ped $element.ped --max-maf 0.01 --kin $element.kinf --test emmaxCMC --out BcA/$element-CMC-0.01
done

理想情况下，这就是我希望脚本执行的操作。

epacts group --groupf CDES-genes.grp --vcf CDES.vcf.gz --ped CDES_MyopV1.ped --max-maf 0.05 --kin CDES.kinf  --test emmaxCMC --out BcA/CDES_MyopV1-CMC-0.05
       epacts group --groupf CDES-genes.grp --vcf CDES.vcf.gz --ped CDES_MyopV1.ped --max-maf 0.03 --kin CDES.kinf  --test emmaxCMC --out BcA/CDES_MyopV1-CMC-0.03 
       epacts group --groupf CDES-genes.grp --vcf CDES.vcf.gz --ped CDES_MyopV1.ped --max-maf 0.02 --kin CDES.kinf --test emmaxCMC --out BcA/CDES_MyopV1-CMC-0.02 
       epacts group --groupf CDES-genes.grp --vcf CDES.vcf.gz --ped CDES_MyopV1.ped --max-maf 0.01 --kin CDES.kinf --test emmaxCMC --out BcA/CDES_MyopV1-CMC-0.01
shell-script bioinformatics
Martin Hope
austin7923
Asked: 2022-04-20 17:02:38 +0800 CST
  • 1

我有一个排序的 ID 和数字（位置）文件。我需要将第二列中的位置分组为一组 500 的间隔。

如果该行的值与上一行相比小于500，则将它们分组到同一组中；而如果该行的值超过 500，则将它们分组到不同的组中。

输入文件：

snp00001    200
snp00002    300
snp00003    400
snp00004    500
snp00005    600
snp00006    900
snp00007    1500
snp00008    1800
snp00009    3000
snp00010    3500
snp00011    4000
snp00012    5000

期望的输出

snp00001 200 Group1
snp00002 300 Group1
snp00003 400 Group1
snp00004 500 Group1
snp00005 600 Group1
snp00006 900 Group1
snp00007 1500 Group2
snp00008 1800 Group2
snp00009 3000 Group3
snp00010 3500 Group3
snp00011 4000 Group4
snp00012 5000 Group5

额外说明：snp00001 到 snp00006 将被归入同一组，因为它们之间的范围 (snp00002 - snp00001) 或 (snp00003 - snp00002) 或 (snp00004 - snp00003) ... 小于 500。

snp00006 和 snp00007 被分到下一组，因为它们之间的范围（snp00007 - snp00006）超过 500。

我试过用awk，但没有成功。

awk -v step=500 -v OFS='\t' '{if(NR==1 || $2+limit){group++} file="Group"group; print file,$0}' input_file
text-processing bioinformatics
Martin Hope
austin7923
Asked: 2022-04-12 00:11:35 +0800 CST
  • 0

我有一个排序的 ID 和数字（位置）文件。我需要将第二列中的位置分组为 500 个间隔，然后拆分为不同的文件。

输入

snp00001    200
snp00002    300
snp00003    400
snp00004    500
snp00005    600
snp00006    900
snp00007    1500
snp00008    1800
snp00009    3000
snp00010    3500
snp00011    4000
snp00012    5000

期望的输出

snp00001    200 Group1
snp00002    300 Group1
snp00003    400 Group1
snp00004    500 Group1
snp00005    600 Group1
snp00006    900 Group2
snp00007    1500    Group3
snp00008    1800    Group3
snp00009    3000    Group4
snp00010    3500    Group4
snp00011    4000    Group5
snp00012    5000    Group6

然后将这些组保存到不同的文件中，分别重命名为Group1、和。Group2Group3Group4

我已经尝试bedtools了一些其他命令，但问题无法解决。

任何帮助将不胜感激。

谢谢！

text-processing bioinformatics
Martin Hope
Callahan McGovern
Asked: 2022-04-06 06:01:55 +0800 CST
  • 1

我想删除文件中每次出现的the>和 word之间的部分。Bacteria

  • 这是输入的示例： 
    >AADV02000003.105686.107093 Bacteria;Cyanobacteria;Cyanobacteriia;Cyanobacteriales;
  • 结果应如下所示： 
    >Bacteria;Cyanobacteria;Cyanobacteriia;Cyanobacteriales;

这是一个 FASTA 文件（仿生信息学中的一种常见文件格式），因此>始终是该行中的第一个非空白字符，并且该行中只有一个这样的字符。

我正在考虑使用sed，但我不确定如何设置命令。感谢您的帮助。

text-processing bioinformatics
Martin Hope
Darker Walker
Asked: 2021-12-14 23:59:23 +0800 CST
  • -1

我有一个包含特定蛋白质的多个序列（氨基酸序列）的 fasta 文件。序列的最后一个字符表示为“*”，它实际上代表终止密码子。我正在尝试使用 MUSCLE 进行多序列比对，但该工具拒绝所有最后带有“*”的序列。

如何删除“*”？

例如我的输入文件是：

>seq1
MSDGFHS*
>Seq2
MSDRFH*

我需要的是：

>seq1
MSDGFHS
>Seq2
MSDRFH
text-processing bioinformatics
Martin Hope
structural_lexa
Asked: 2021-11-11 14:15:01 +0800 CST
  • 2

我有大约 4500 行氨基酸变异，如下所示：

S1437T
H1266Y
T2662A
E1397A
E626K
S1538T
E3021K

简而言之，数字两侧的字母是氨基酸残基，数字代表残基位置。我只想检索那些在 2400 到 3100 范围内的变化。

我尝试使用grep，但没有那么成功。另外，我知道这awk对于这种操作可能会更好，但在awk. 任何帮助，将不胜感激。

text-processing bioinformatics
Martin Hope
Paulo Sergio Schlogl
Asked: 2021-10-22 03:38:34 +0800 CST
  • 1

我发现了这个有趣的命令：

grep -v '^>' test.fasta | tr -d '\n' | sed -e 's/\(.\)/\1\n/g' | sort | uniq -c | sort -rn

我对它的含义有所了解（它计算文本文件中的字母），但我的问题是关于这个：

sed -e 's/\(.\)/\1\n/g'

我知道它由三个替代命令组成。一种是替换换行符（\n），一个匹配除换行符（\(.\)）之外的任何字符，但我迷路了/\1\？

sed bioinformatics
Martin Hope
moth
Asked: 2021-10-21 00:42:35 +0800 CST
  • 2

我有一个看起来像这样的文件：

chr1    3143567 3143568 .3-2704 1.000000|ENSMUSG00000102693.2
chr1    3143599 3143600 .3-2705 1.000000|ENSMUSG00000102693.2
chr1    3143631 3143632 .3-2706 1.000000|ENSMUSG00000102693.2
chr1    3143663 3143664 .3-2707 1.000000|ENSMUSG00000102693.2
chr1    3143695 3143696 .3-2708 1.000000|ENSMUSG00000102693.2
chr1    3143727 3143728 .3-2709 1.000000|ENSMUSG00000102693.2

我正在编写 2 个 sed 表达式来过滤|第一个之前的所有内容，并使用结果文件我丢弃之后的所有内容，.如下所示：

sed -n -e 's/^.*|//p' original_file.txt > first_result.txt

sed -n -e 's/\..*//p' first_result.txt > final_result.txt

我怎样才能将所有这些写在一行中？

最终目标是捕获ENSMUSG00000102693

sed bioinformatics
Martin Hope
Ricardo Guerreiro
Asked: 2020-01-25 07:56:27 +0800 CST
  • 6

所以我有这个fasta（生物学）文件，看起来像这样：

>m64093_191209_130050/133911/ccs_64
TTCAGGCTGTGTTCCATTTGATTTAAAATCAAATAATTTCATTCGCGTCAGAACACCTGGTTTCACGACC
ATAAATAATTTACCAGTGAATCGAGGCTCAATTATAGATCCTCGGACGCGAGTTCTCGGTTGACGAGTGG
GATTCGAATTATTTTTCACCGAAAATTTTAGTCGACGAGTTCAGATAAATTTGTTCGGGATAAAATCATC
TGAGTAGGTCGGGCTTCTGAATTTCGTATTCTTGCGAGCAATGAATTTTAAATAATCATCGGACATACCA
ATTTTTGGAACAATAATGTTCCGAACATCCCGAAAATATAGGAAGAGCCCGGATAGATAAAAATAAACAC

每行最长为 70 个字符。通常，如果我想将其格式化为最多 50 个字符，我使用：

fold -50 input.fasta > output.fasta # 也试过 -b 和 -w args

但不知何故，这是行不通的。该文件看起来与我见过的许多其他文件完全相同。输出现在如下所示：

>m64093_191209_130050/133911/ccs_64
TTCAGGCTGTGTTCCATTTGATTTAAAATCAAATAATTTCATTCGCGTCA
GAACACCTGGTTTCACGACC
ATAAATAATTTACCAGTGAATCGAGGCTCAATTATAGATCCTCGGACGCG
AGTTCTCGGTTGACGAGTGG
GATTCGAATTATTTTTCACCGAAAATTTTAGTCGACGAGTTCAGATAAAT
TTGTTCGGGATAAAATCATC
TGAGTAGGTCGGGCTTCTGAATTTCGTATTCTTGCGAGCAATGAATTTTA
AATAATCATCGGACATACCA
ATTTTTGGAACAATAATGTTCCGAACATCCCGAAAATATAGGAAGAGCCC

它剪切了突出的 20 个字符并将它们正确放置在下面，但是它没有加入下一行并将其剪切到最多 50 个字符上。

我回到以前创建的 fasta 文件， fold 命令仍然正常工作。如果我复制新文件的一部分并将其粘贴到另一个文件中，问题仍然存在。

我认为可能存在我不知道的编码问题。任何人都可以帮忙吗？

干杯，

编辑：很好的答案，谢谢！

bioinformatics fold