如何在R语言中拆分鸟枪法基因组FASTA文件的编码与非编码区

Alright, let's walk through how to split your shotgun genome sequence (205,000 bp, accession NZ_LRPF01000001) into coding (CDS) and non-coding regions using R. I'll break this down into straightforward steps you can follow easily:

We'll use Bioconductor's Biostrings for sequence manipulation and rentrez to pull the sequence directly from NCBI. Install them if you haven't already:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("Biostrings", "rentrez"))

Load the packages and retrieve your sequence using the accession ID (no need to navigate to external URLs manually):

library(Biostrings)
library(rentrez)

# Fetch the FASTA record from NCBI
seq_record <- entrez_fetch(db = "nuccore", id = "NZ_LRPF01000001", rettype = "fasta")
# Convert to a DNAStringSet object for easy manipulation
genome_seq <- readDNAStringSet(textConnection(seq_record))

Note: If you already have the FASTA file saved locally, replace the above with genome_seq <- readDNAStringSet("your_local_sequence.fasta").

Based on your data, we'll formalize the CDS intervals (remember these are 1-based positions, which matches standard genomic formatting):

# Create a data frame with CDS start and end positions
cds_regions <- data.frame(
    start = c(343, 937),
    end = c(780, 1866)
)

Use extractAt() to pull out the CDS regions from the full genome sequence:

# Convert CDS positions to an IRanges object
cds_ranges <- IRanges(start = cds_regions$start, end = cds_regions$end)
# Extract CDS sequences
cds_seqs <- extractAt(genome_seq, cds_ranges)
# Rename sequences for clarity
names(cds_seqs) <- paste0("CDS_", 1:nrow(cds_regions))

To get non-coding regions, we'll subtract the CDS ranges from the full genome range. This automatically captures all the non-CDS intervals you mentioned (1-342, 781-936, 1867-205000):

# Define the full range of the genome
full_genome_range <- IRanges(start = 1, end = width(genome_seq))
# Get non-coding ranges by removing CDS regions from the full genome
non_cds_ranges <- setdiff(full_genome_range, cds_ranges)
# Extract non-coding sequences
non_cds_seqs <- extractAt(genome_seq, non_cds_ranges)
# Rename non-coding sequences
names(non_cds_seqs) <- paste0("NonCDS_", 1:length(non_cds_ranges))

If you want to export the sequences to FASTA files for downstream analysis:

# Save coding sequences
writeXStringSet(cds_seqs, "coding_sequences.fasta")
# Save non-coding sequences
writeXStringSet(non_cds_seqs, "non_coding_sequences.fasta")

• Double-check your CDS start/end positions to ensure they're 1-based (this is standard for genomic sequences from public databases).
• If you need to add more CDS regions later, just extend the cds_regions data frame with new start/end pairs.

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