统计file.txt中唯一DNA序列出现次数及低内存实现方案咨询

Hey there! Let's break down how to count your DNA sequences efficiently, while keeping memory usage in check—since you’ve got tons of these files to work with later. Here’s a step-by-step guide tailored to your needs:

1. Generate Your Desired Outputs #

First, let’s cover both the single-column matrix and CSV output options, starting with methods that balance speed and memory efficiency.

Option 1: Single-Column Matrix (Row Names = Unique Sequences) #

We’ll use scan (more memory-efficient than readLines for large text files) to load sequences, then count and convert to a matrix:

# Load sequences without extra overhead
dna_seqs <- scan("file.txt", what = character(), quiet = TRUE)

# Count occurrences and convert to a matrix
count_table <- table(dna_seqs)
count_matrix <- as.matrix(count_table)

# Clean up labels
colnames(count_matrix) <- "Abundance"
rownames(count_matrix) <- names(count_table)

Option 2: CSV Output (Sequence Headers + Counts) #

For a CSV with unique sequences as the first row and counts as the second, we’ll use data.table (fast, low-memory file handling):

library(data.table)

# Read the file quickly with minimal memory usage
dt <- fread("file.txt", header = FALSE, col.names = "Sequence")

# Count sequences per unique entry
count_dt <- dt[, .(Abundance = .N), by = Sequence]

# Reshape to wide format (matches your desired CSV structure)
wide_dt <- dcast(count_dt, . ~ Sequence, value.var = "Abundance")

# Write to CSV (no row names to keep the format clean)
fwrite(wide_dt, "dna_counts.csv", row.names = FALSE)

2. Memory Optimization for Large Numbers of Files #

Since you’ll need to merge results across many files, minimizing memory usage during individual file processing is key. Here are the most impactful strategies:

a. Process Files in Chunks (Avoid Loading All Data at Once) #

Instead of loading the entire file into memory, read it in chunks, count sequences per chunk, and accumulate results. This keeps memory usage low even for massive files:

# Open a connection to the file
con <- file("file.txt", "r")
total_counts <- integer(0)

# Read 10,000 lines at a time (adjust chunk size based on your system)
while (length(seqs_chunk <- readLines(con, n = 10000)) > 0) {
  chunk_counts <- table(seqs_chunk)
  
  # Merge chunk counts into total counts
  for (seq_name in names(chunk_counts)) {
    total_counts[seq_name] <- ifelse(is.na(total_counts[seq_name]), 
                                     chunk_counts[seq_name], 
                                     total_counts[seq_name] + chunk_counts[seq_name])
  }
}

# Close the file connection
close(con)

# Convert to your desired matrix format
count_matrix <- as.matrix(total_counts)
colnames(count_matrix) <- "Abundance"

b. Use Efficient Data Structures for Counting #

Hash tables (from the hash package) are faster and more memory-efficient than base R vectors for dynamic counting:

library(hash)

con <- file("file.txt", "r")
count_hash <- hash()

while (length(seqs_chunk <- readLines(con, n = 10000)) > 0) {
  chunk_counts <- table(seqs_chunk)
  
  for (seq_name in names(chunk_counts)) {
    if (has.key(seq_name, count_hash)) {
      count_hash[[seq_name]] <- count_hash[[seq_name]] + chunk_counts[seq_name]
    } else {
      count_hash[[seq_name]] <- chunk_counts[seq_name]
    }
  }
}

close(con)

# Convert hash to matrix
count_matrix <- as.matrix(unlist(count_hash))
rownames(count_matrix) <- names(count_hash)
colnames(count_matrix) <- "Abundance"

c. Merge Results Without Re-loading Raw Data #

When combining results across multiple files, only store the count data (not the raw sequences) to save memory. Here’s how to merge with data.table:

library(data.table)

# List all your DNA sequence files
file_list <- list.files(path = "your_directory", pattern = "*.txt", full.names = TRUE)

# Initialize an empty table to hold combined counts
all_counts <- data.table(Sequence = character(), Abundance = integer())

for (file in file_list) {
  # Read and count sequences in the current file
  dt <- fread(file, header = FALSE, col.names = "Sequence")
  chunk_counts <- dt[, .(Abundance = .N), by = Sequence]
  
  # Merge with total counts, filling missing values with 0
  all_counts <- merge(all_counts, chunk_counts, by = "Sequence", all = TRUE, suffixes = c("", "_new"))
  all_counts[, Abundance := ifelse(is.na(Abundance), 0, Abundance) + ifelse(is.na(Abundance_new), 0, Abundance_new)]
  all_counts[, Abundance_new := NULL]
}

# Convert to matrix if needed
final_matrix <- as.matrix(all_counts$Abundance)
rownames(final_matrix) <- all_counts$Sequence
colnames(final_matrix) <- "Total_Abundance"

d. Encode Sequences as Factors #

If your DNA sequences are of fixed length, converting them to factors reduces memory usage significantly (factors store unique values once, with pointers to each occurrence):

dna_seqs <- scan("file.txt", what = character(), quiet = TRUE)
dna_seqs <- factor(dna_seqs)  # This cuts memory usage for large datasets
count_matrix <- as.matrix(table(dna_seqs))
colnames(count_matrix) <- "Abundance"

3. Bonus: Ultra-Large Files? Try the ff Package #

For files too big to handle even in chunks, the ff package stores data on disk instead of memory. It’s a bit more complex, but perfect for extreme cases.

内容的提问来源于stack exchange，提问作者colin