标签: bioinformatics
FASTA算法解释
我试图理解FASTA算法在数据库中搜索查询序列的类似序列的基本步骤.这些是算法的步骤:
- 识别I和J之间的常见k字
- 使用k字匹配对角线进行评分,确定10个最佳对角线
- 使用替换分数矩阵重新构造初始区域
- 使用间隙加入初始区域,惩罚差距
- 执行动态编程以查找最终对齐
我对使用PAM250评分矩阵的第3步和第4步以及如何"加入使用差距"感到困惑.
有人可以"尽可能具体地"为我解释这两个步骤.谢谢
推荐指数
解决办法
查看次数
使用元组填充列表
我只是在模仿(孟德尔的第一遗传定律).
在我可以让生物交配并分析结果之前,必须生成群体,即,必须在不解包它们的情况下用不同数量的三种不同类型的元组填充列表.
在尝试熟悉itertools时(我稍后需要在配合部分中进行组合),我想出了以下解决方案:
import itertools
k = 2
m = 3
n = 4
hd = ('A', 'A') # homozygous dominant
het = ('A', 'a') # heterozygous
hr = ('a', 'a') # homozygous recessive
fhd = itertools.repeat(hd, k)
fhet = itertools.repeat(het, m)
fhr = itertools.repeat(hr, n)
population = [x for x in fhd] + [x for x in fhet] + [x for x in fhr]
这将导致:
[('A', 'A'), ('A', 'A'), ('A', 'a'), ('A', 'a'), ('A', 'a'), …推荐指数
解决办法
查看次数
使用Brainarray定制CDF对GSE31312的RMA标准化中的所有NaN进行标准化
我正在尝试使用来自Brainarray的定制基因水平注释CDF(芯片定义文件)文件,使关于弥漫性大B细胞淋巴瘤的特定基因表达数据集正常化.
不幸的是,RMA归一化表达式矩阵是所有NaNs,我不明白为什么.
数据集(GSE31312)可在GEO网站免费获取,并使用Affymetrix HG-U133 Plus 2.0阵列平台.我正在使用该affy包来执行RMA规范化.
由于问题是特定于数据集的,因此遗憾的是,以下用于重现问题的R代码非常麻烦(2 GB下载,8.8 GB解压缩).
设置工作目录.
setwd("~/Desktop/GEO")
加载所需的包.取消注释以安装软件包.
#source("http://bioconductor.org/biocLite.R")
#biocLite(pkgs = c("GEOquery", "affy", "AnnotationDbi", "R.utils"))
library("GEOquery") # To automatically download the data
library("affy")
library("R.utils") # For file handling
将阵列数据下载到工作目录.
files <- getGEOSuppFiles("GSE31312")
解压缩名为CEL的目录中的数据
#Sys.setenv(TAR = '/usr/bin/tar') # For (some) OS X uncommment this line
untar(tarfile = "GSE31312/GSE31312_RAW.tar", exdir = "CEL")
解压缩所有.gz文件
gz.files <- list.files("CEL", pattern = "\\.gz$",
ignore.case = TRUE, full.names = TRUE)
for (file in …推荐指数
解决办法
查看次数
从仅包含0或仅包含0的数据框中删除行
我正在尝试在R中创建一个函数,它允许我根据行是否包含一个零的单个列来过滤我的数据集.此外,有时我只想删除所有列中为零的行.
而且,这是它变得有趣的地方; 并非所有列都包含数字,列数可能会有所不同.
我试图将我的一些数据粘贴到我想要获得的结果中.
unfiltered:
ID GeneName DU145small DU145total PC3small PC3total
1 MIR22HG 33221.5 1224.55 2156.43 573.315
2 MIRLET7E 87566.1 7737.99 25039.3 16415.6
3 MIR612 0 0 530.068 0
4 MIR218-1 0 0 1166.88 701.253
5 MIR181B2 70723.2 3958.01 6209.85 1399.34
6 MIR218-2 0 0 0 0
7 MIR10B 787.516 330.556 0 20336.4
8 MIR3176 0 0 0 0
any rows with containing a zero removed:
ID GeneName DU145small DU145total PC3small PC3total
1 MIR22HG 33221.5 1224.55 2156.43 573.315
2 MIRLET7E …推荐指数
解决办法
查看次数
在滑动窗口中找到k-mers
我正在尝试解决这个生物信息学问题:https://stepic.org/lesson/An-Explosion-of-Hidden-Messages-4/step/1?course= Bioinformatics-Algorithms-2 &unit=8
具体问题在上面链接的第5个窗口,问题是: 大肠杆菌基因组中有多少个不同的9聚体形成(500,3)个团块？(换句话说,不要多次计算9-mer.)
我的代码如下.这是错误的,我想要解释为什么,以及如何改进它(显然O效率很糟糕,但几天前我开始编写Python ...)非常感谢!
genome = '' #insert e. Coli genome here
k = 4 #length of k-mer
L = 50 #size of sliding window
t = 3 #k-mer appears t times
counter = 0
Count = []
for i in range(0,len(genome)-L): #slide window down the genome
pattern = genome[i:i+k] #given this k-mer
for j in range(i,i+L): #calculate k-mer frequency in window of len(L)
if genome[j:j+k] == pattern:
counter = counter + 1
Count.append(counter) …推荐指数
解决办法
查看次数
prcomp和ggbiplot:无效的'rot'值
我正在尝试使用R对我的数据进行PCA分析,我找到了这个很好的指南,使用prcomp和ggbiplot.我的数据是两种样本类型,每种类型有三个生物重复(即6行)和大约20000个基因(即变量).首先,使用指南中描述的代码获取PCA模型不起作用:
>pca=prcomp(data,center=T,scale.=T)
Error in prcomp.default(data, center = T, scale. = T) :
cannot rescale a constant/zero column to unit variance
但是,如果我删除该scale. = T部件,它工作正常,我得到一个模型.这是为什么,这是下面错误的原因？
> summary(pca)
Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation 4662.8657 3570.7164 2717.8351 1419.3137 819.15844
Proportion of Variance 0.4879 0.2861 0.1658 0.0452 0.01506
Cumulative Proportion 0.4879 0.7740 0.9397 0.9849 1.00000
其次,绘制PCA.即使只是使用基本代码,我得到一个错误和一个空的情节图像:
> ggbiplot(pca)
Error: invalid 'rot' value
这意味着什么,我该如何解决？是否与制作PCA的(非)规模有关,还是有所不同？我认为它必须与我的数据有关,因为如果我使用标准示例代码(下面),我会得到一个非常好的PCA图.
> data(wine)
> wine.pca=prcomp(wine,scale.=T)
> print(ggbiplot(wine.pca, obs.scale = 1, …推荐指数
解决办法
查看次数
使用Gibbs采样器搜索主题
我是编程和生物信息学的初学者.所以,我很感激你的理解.我尝试使用Gibbs采样开发一个用于motif搜索的python脚本,如Coursera类"在DNA中查找隐藏的消息"中所述.课程中提供的伪代码为:
GIBBSSAMPLER(Dna, k, t, N)
randomly select k-mers Motifs = (Motif1, …, Motift) in each string
from Dna
BestMotifs ? Motifs
for j ? 1 to N
i ? Random(t)
Profile ? profile matrix constructed from all strings in Motifs
except for Motifi
Motifi ? Profile-randomly generated k-mer in the i-th sequence
if Score(Motifs) < Score(BestMotifs)
BestMotifs ? Motifs
return BestMotifs
问题描述:
代码挑战:实施GIBBSSAMPLER.
输入:整数k,t和N,后跟Dna的字符串集合.输出:由20个随机启动运行GIBBSSAMPLER(Dna,k,t,N)产生的字符串BestMotifs.记得使用伪版!
样本输入:
8 5 100
CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG
TAGTACCGAGACCGAAAGAAGTATACAGGCGT
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC
AATCCACCAGCTCCACGTGCAATGTTGGCCTA
样本输出:
TCTCGGGG
CCAAGGTG
TACAGGCG …推荐指数
解决办法
查看次数
来自gff3文件的scikit-bio extract基因组特征
scikit-bio是否有可能从基因组fasta文件中提取存储在gff3格式文件中的基因组特征？
例:
genome.fasta
>sequence1
ATGGAGAGAGAGAGAGAGAGGGGGCAGCATACGCATCGACATACGACATACATCAGATACGACATACTACTACTATGA
annotation.gff3
#gff-version 3
sequence1 source gene 1 78 . + . ID=gene1
sequence1 source mRNA 1 78 . + . ID=transcript1;parent=gene1
sequence1 source CDS 1 6 . + 0 ID=CDS1;parent=transcript1
sequence1 source CDS 73 78 . + 0 ID=CDS2;parent=transcript1
mRNA特征(转录物1)的所需序列将是两个子CDS特征的连接.所以在这种情况下,这将是'ATGGAGCTATGA'.
推荐指数
解决办法
查看次数
如何着色绘图子区域并使用ggrepel标记数据点的子集？
我制作了这个火山阴谋,并希望改进它如下:
用蓝色数据点完全遮蔽该区域:使用我当前的代码,我无法将阴影扩展到您所看到的范围之外.我希望它一直到剧情区域限制.
geom_text让我来标记数据点的子集,但做起来ggrepel应该添加有标签从而提高标签的清晰度连接数据点的线.如何重用现有geom_text代码ggrepel来实现这一目标？
这是我的代码:
ggplot(vol.new, aes(x = log2.fold.change, y = X.NAME., fill = Color)) + # Define data frame to be used for plotting; define data for x and y axes; crate a scatterplot object.
geom_point(size = 2, shape = 21, colour = "black") + # Define data point style.
ggtitle(main.title, subtitle = "Just a little subtitle") + # Define title and subtitle.
labs(x = x.lab, y = y.lab) + # Define labels …推荐指数
解决办法
查看次数
通过对列进行分组将数据帧拆分为多个数据帧
我有一个表达数据的数据框，其中基因是行，列是样本。我还有一个数据框，其中包含表达式数据框中每个样本的元数据。实际上，我的 expr 数据框有 30,000 多行和 100 多列。然而，下面是一个数据较小的示例。
expr <- data.frame(sample1 = c(1,2,2,0,0),
sample2 = c(5,2,4,4,0),
sample3 = c(1,2,1,0,1),
sample4 = c(6,5,6,6,7),
sample5 = c(0,0,0,1,1))
rownames(expr) <- paste0("gene",1:5)
meta <- data.frame(sample = paste0("sample",1:5),
treatment = c("control","control",
"treatment1",
"treatment2", "treatment2"))
我想找到每次治疗中每个基因的平均值。从我看到的 split() 或 group_by() 示例中，人们根据 data.frame 中已存在的列进行分组。但是，我有一个单独的数据框（元），用于对另一个数据框（expr）中的列进行分组。
我希望我的输出是一个数据框，其中基因作为行，治疗作为列，值作为平均值。
# control treatment1 treatment2
# gene1 mean mean mean
# gene2 mean mean mean
推荐指数
解决办法
查看次数