我有以下任务:从30个字符长的模式序列开始(它实际上是DNA序列,以免称它为P30)我需要在文本文件中找到所有行开始(^ agacatacag ...)并使用精确的P30,然后使用29最后一个字符是30,28和最多10个字符.我需要做的就是删除模式的第一个字符并继续搜索.为简单起见,我目前需要精确匹配,但允许1个不匹配更长时间(20-30个字符长的模式)会更好.
我目前的一个相当慢的解决方案是创建一个每行有一个截断模式的shell文件,然后grep [1].这意味着我正在阅读巨大的,少量的GB文本文件20x,这可能需要一天+.
我可以切换到python,创建一个包含所有必需模式的列表/元组,然后只读取一次文件,每个序列循环20次,使用pypy加快速度.
[1]因为它是DNA序列,并且要搜索的输入是FASTQ格式,我使用fqgrep:https://github.com/indraniel/fqgrep with tre library:https://github.com/laurikari/tre/
edit_1 更改示例(缩短模式).只显示前几步/更短的模式:
^abcde
^bcde
^cde
或者如果你喜欢它作为DNA:
^GATACCA
^ATACCA
^TACCA
edit_2 简单的grep并没有真正削减它.我需要对每行只有第2行匹配的4行FASTQ格式进行后处理.如果我不使用fqgrep,那么我必须:
读取输入的4行
- 检查第2行(序列)是否以20种模式中的任何一种开始(P30-P10)
- 如果我得到了匹配,我需要删除第2行的第N个字符和# 4,其中N代表匹配模式的长度 - 打印输出/写入文件行#1- $ 4无匹配无效
对于内部解决方案,我可以尝试使用GNU并行分割输入文件,例如4M的谎言块,并以这种方式加快速度.但是,如果我想让其他人使用每个新软件,我要求最终用户安装广告,这是一个额外的复杂程度.
**编辑3**来自Vyctor的正则表达式和匹配行的简单示例:
starting P30 regex
^agacatacagagacatacagagacatacag
matching sequence:
^agacatacagagacatacagagacatacagGAGGACCA
P29:
^gacatacagagacatacagagacatacag
matching sequence:
^gacatacagagacatacagagacatacagGACCACCA
P28:
^acatacagagacatacagagacatacag
matching sequence:
^acatacagagacatacagagacatacagGATTACCA
我请从左侧的字符/ DNA碱基(或在DNA讲5'端),因为这是这些序列是由真实的酶降解的方式.一旦找到正则表达式序列本身就没有意义.所需的输出是正则表达式之后的读取序列.在上面的例子中,它在UPERCASE中,然后可以在下一步中映射到基因组.应该强调的是,除了这个玩具示例之外,我正在变得更长,在正则表达式之后的先验未知和变化的序列.在现实世界中,我不必处理DNA的大小写字符(一切都是大写的),但我可能会在我搜索模式的序列中遇到Ns(=未知的DNA碱基).这些可以在第一次近似中被忽略,但是对于更敏感的算法版本,可能应该以简单的不匹配来处理.在理想情况下,不考虑给定位置的简单不匹配,而是计算更复杂的罚分,同时考虑存储在FASTQ格式中存储的每4行长序列记录的第4行中的DNA序列质量值:http://en.wikipedia .ORG /维基/ FASTQ_format#质量
但这种方式更为复杂,到目前为止,"只采用与正则表达式完美匹配的读取"方法已经足够好,并使后续步骤更容易分析.
我得到了以下数据框,df,(此处显示的片段):
H2475 H2481 H2669 H2843 H2872 H2873 H2881 H2909
E1 94.470 26.481 15.120 18.490 16.189 11.422 14.886 0.512
E2 1.016 0.363 0.509 1.190 1.855 0.958 0.771 0.815
E3 9.671 0.637 0.571 0.447 0.116 0.452 0.403 0.003
E4 3.448 2.826 2.183 2.607 4.288 2.526 2.820 3.523
E5 2.548 1.916 1.126 1.553 1.089 1.228 0.887 1.065
我想要做的是在删除两个极值后计算每行的平均值.对于整行,我使用了plyr:
library(plyr)
df.my_means <- adply(df, 1, transform, my_means = mean(as.matrix(df[i,]) ) )
创建一些临时数据帧/矩阵也应该可以用最小值和最大值替换为NA,但作为初学者,我无法做到.
非常感谢你的帮助
编辑1
我显然没有意识到这意味着有一个修剪选项.我想有一个解决方案,而不是意味着我可以插入任何其他功能.即:
library(plyr)
library(e1071)
df.my_means <- adply(df, 1, …我得到了以下数据框(这里简化):
H2475 H2481 H2669 H2843 H2872 H2873 H2881 H2909
E1 24.470 26.481 15.120 18.490 16.189 11.422 14.886 18.512
E2 1.016 0.363 0.509 1.190 1.855 0.958 0.771 0.815
E3 0.671 0.637 0.571 0.447 0.116 0.452 0.403 0.753
E4 3.448 2.826 2.183 2.607 4.288 2.526 2.820 3.523
E5 2.548 1.916 1.126 1.553 1.089 1.228 0.887 1.065
其中E1 - E5是行号.我想为每行的值创建一个boxplot(甚至更好的小提琴图).我想忽略列ID,并将行ID作为"因子".
像这样的回答: 如何生成boxplot
但使用E1 - E5作为值.正如你可能猜到我是新手融化/演员/重塑.
非常感谢您的帮助