从排序数据中解析信息 - python

我有一个txt文件，该文件是经过转换的fasta文件，只是具有我要分析的特定区域。看起来像这样

CTGGCCGCGCTGACTCCTCTCGCT

CTCGCAGCACTGACTCCTCTTGCG

CTAGCCGCTCTGACTCCGCTAGCG

CTCGCTGCCCTCACACCTCTTGCA

CTCGCAGCACTGACTCCTCTTGCG

CTCGCAGCACTAACACCCCTAGCT

CTCGCTGCTCTGACTCCTCTCGCC

CTGGCCGCGCTGACTCCTCTCGCT

我目前正在使用excel对每个位置的核苷酸多样性进行一些计算。有些文件读取了200,000次，因此这使excel文件变得笨拙。我认为必须有一种使用python或R的简便方法。

基本上，我想将.txt文件与序列列表一起使用，并使用此方程–p（log2（p））测量每个位置的核苷酸多样性。有谁知道如何用excel以外的方法来做到这一点？

提前非常感谢您的帮助。

参考方案

如果可以使用fasta文件工作，那可能会更好，因为
专为使用该格式而设计的软件包。

在这里，我使用包seqinr在R中提供了一个解决方案，
dplyr（tidyverse的一部分）用于处理数据。

如果这是您的fasta文件（根据您的序列）：

>seq1
CTGGCCGCGCTGACTCCTCTCGCT
>seq2
CTCGCAGCACTGACTCCTCTTGCG
>seq3
CTAGCCGCTCTGACTCCGCTAGCG
>seq4
CTCGCTGCCCTCACACCTCTTGCA
>seq5
CTCGCAGCACTGACTCCTCTTGCG
>seq6
CTCGCAGCACTAACACCCCTAGCT
>seq7
CTCGCTGCTCTGACTCCTCTCGCC
>seq8
CTGGCCGCGCTGACTCCTCTCGCT

您可以使用seqinr包将其读入R：

# Load the packages
library(tidyverse) # I use this package for manipulating data.frames later on
library(seqinr)

# Read the fasta file - use the path relevant for you
seqs <- read.fasta("~/path/to/your/file/example_fasta.fa")

这将返回一个list对象，其中包含的元素数与
文件中的序列。

对于您的特定问题-计算每个职位的多样性指标-
我们可以使用seqinr包中的两个有用的函数：

getFrag()子集化序列
count()计算每个核苷酸的频率

例如，如果我们想要第一个位置的核苷酸频率
我们的序列，我们可以做：

# Get position 1
pos1 <- getFrag(seqs, begin = 1, end = 1)

# Calculate frequency of each nucleotide
count(pos1, wordsize = 1, freq = TRUE)

a c g t 
0 1 0 0 

向我们表明，第一个位置仅包含“ C”。

下面是一种编程方式“循环”所有位置并执行
我们可能对以下计算感兴趣：

# Obtain fragment lenghts - assuming all sequences are the same length!
l <- length(seqs[[1]])

# Use the `lapply` function to estimate frequency for each position
p <- lapply(1:l, function(i, seqs){
  # Obtain the nucleotide for the current position
  pos_seq <- getFrag(seqs, i, i)

  # Get the frequency of each nucleotide
  pos_freq <- count(pos_seq, 1, freq = TRUE)

  # Convert to data.frame, rename variables more sensibly
  ## and add information about the nucleotide position
  pos_freq <- pos_freq %>% 
    as.data.frame() %>%
    rename(nuc = Var1, freq = Freq) %>% 
    mutate(pos = i)
}, seqs = seqs)

# The output of the above is a list.
## We now bind all tables to a single data.frame
## Remove nucleotides with zero frequency
## And estimate entropy and expected heterozygosity for each position
diversity <- p %>% 
  bind_rows() %>% 
  filter(freq > 0) %>% 
  group_by(pos) %>% 
  summarise(shannon_entropy = -sum(freq * log2(freq)),
            het = 1 - sum(freq^2), 
            n_nuc = n())

这些计算的输出现在看起来像这样：

head(diversity)

# A tibble: 6 x 4
    pos shannon_entropy     het n_nuc
  <int>           <dbl>   <dbl> <int>
1     1        0.000000 0.00000     1
2     2        0.000000 0.00000     1
3     3        1.298795 0.53125     3
4     4        0.000000 0.00000     1
5     5        0.000000 0.00000     1
6     6        1.561278 0.65625     3

这是它的更直观视图（使用ggplot2，它也是tidyverse包的一部分）：

ggplot(diversity, aes(pos, shannon_entropy)) + 
  geom_line() +
  geom_point(aes(colour = factor(n_nuc))) +
  labs(x = "Position (bp)", y = "Shannon Entropy", 
       colour = "Number of\nnucleotides")

更新：

要将其应用于多个fasta文件，这是一种可能性
（我没有测试此代码，但是类似的东西应该可以工作）：

# Find all the fasta files of interest
## use a pattern that matches the file extension of your files
fasta_files <- list.files("~/path/to/your/fasta/directory", 
                          pattern = ".fa", full.names = TRUE)

# Use lapply to apply the code above to each file
my_diversities <- lapply(fasta_files, function(f){
  # Read the fasta file
  seqs <- read.fasta(f)

  # Obtain fragment lenghts - assuming all sequences are the same length!
  l <- length(seqs[[1]])

  # .... ETC - Copy the code above until ....
  diversity <- p %>% 
    bind_rows() %>% 
    filter(freq > 0) %>% 
    group_by(pos) %>% 
    summarise(shannon_entropy = -sum(freq * log2(freq)),
              het = 1 - sum(freq^2), 
              n_nuc = n())
})

# The output is a list of tables. 
## You can then bind them together, 
## ensuring the name of the file is added as a new column "file_name"

names(my_diversities) <- basename(fasta_files) # name the list elements
my_diversities <- bind_rows(my_diversities, .id = "file_name") # bind tables

这将为您提供每个文件的多样性表。然后，您可以使用ggplot2对其进行可视化，类似于我上面所做的，但也许可以使用facets将每个文件的多样性分为不同的面板。

Python Pandas导出数据 - python

我正在使用python pandas处理一些数据。我已使用以下代码将数据导出到excel文件。writer = pd.ExcelWriter('Data.xlsx'); wrong_data.to_excel(writer,"Names which are wrong", index = False); writer.…

Python:检查是否存在维基百科文章 - python

我试图弄清楚如何检查Wikipedia文章是否存在。例如，https://en.wikipedia.org/wiki/Food 存在，但是https://en.wikipedia.org/wiki/Fod 不会，页面只是说：“维基百科没有此名称的文章。”谢谢！ 参考方案 >>> import urllib >>> prin…

Python GPU资源利用 - python

我有一个Python脚本在某些深度学习模型上运行推理。有什么办法可以找出GPU资源的利用率水平？例如，使用着色器，float16乘法器等。我似乎在网上找不到太多有关这些GPU资源的文档。谢谢！ 参考方案 您可以尝试在像Renderdoc这样的GPU分析器中运行pyxthon应用程序。它将分析您的跑步情况。您将能够获得有关已使用资源，已用缓冲区，不同渲染状态上…

Python pytz时区函数返回的时区为9分钟 - python

由于某些原因，我无法从以下代码中找出原因：>>> from pytz import timezone >>> timezone('America/Chicago') 我得到：<DstTzInfo 'America/Chicago' LMT-1 day, 18:09:00 STD…

Python uuid4，如何限制唯一字符的长度 - python

在Python中，我正在使用uuid4（）方法创建唯一的字符集。但是我找不到将其限制为10或8个字符的方法。有什么办法吗？uuid4()ffc69c1b-9d87-4c19-8dac-c09ca857e3fc谢谢。 参考方案 尝试：x = uuid4() str(x)[:8] 输出："ffc69c1b" Is there a way to…