Hello,

Thanks.

The answers to your questions depend on what type of analysis you are trying to perform. Essentially, dNdScv uses a reference genome and a table of mutations with positions mapped to the reference. For most applications in next generation sequencing, users map reads to the reference genome and call mutations on them, in which case the idea of a reference genome and a table of mutations is commonplace.

For SARS-CoV-2, you typically have a fasta file with multiple genomes. The way you run dNdScv depends on the questions that you want to address. For example, I wanted to calculate dN/dS on the mutations that have occurred between the Wuhan reference and Omicron. I had two Omicron genomes and the MN908947.3 reference (Wuhan). I used MAFFT to align the two Omicron genomes against the reference using a command like this:

mafft.bat --auto --thread -1 --addfragments omicron.fa MN908947.3.reference.fasta > MSA_with_reference3.fa

This generates a multiple sequence alignment with the 3 genomes, with MN908947.3 at the top.

MSA_with_reference3.fa.txt

I then ran a custom R script to generate a table of mutations for dNdScv listing all of the differences between the two Omicron genomes and the reference.

fa = seqinr::read.fasta("MSA_with_reference3.fa.txt")

mutations = NULL
for (j in 1:length(fa)) {
  r = toupper(fa[[1]])
  s = toupper(fa[[j]])
  mms = which(r != s & s != "N")
  if (any(mms)) {
    refpos = cumsum(r != "-") # This is unnecessary here, but it serves to handle insertions (shows as gaps in the reference) in other MSAs
    m = data.frame(sampleID = names(fa)[j], chr = "MN908947.3", pos = refpos[mms], ref = r[mms], mut = s[mms], stringsAsFactors = F)
    mutations = rbind(mutations, m)
  }
}

# Collapsing the mutation table into unique mutations

mutations$sampleID = "Unique"
mutations = unique(mutations)

# Grouping consecutive mutations into MNVs

annotate_runs = function(m) {
  m = m[order(m$sampleID,m$pos), ]
  runs = rle(m$pos-(1:nrow(m)))
  rmrow = rep(F,nrow(m))
  runstarts = cumsum(runs$length)-runs$length+1
  for (h in 1:length(runs$length)) {
    if (runs$length[h]>1) { # Adjacent mutations
      mutcluster = runstarts[h]:(runstarts[h]+runs$lengths[h]-1)
      rmrow[mutcluster[-1]] = T # Removing all the affected rows except the first one, where we will annotate the MNV
      m[mutcluster[1],"ref"] = paste(m[mutcluster,"ref"],collapse="")
      m[mutcluster[1],"mut"] = paste(setdiff(m[mutcluster,"mut"],"-"),collapse="")
    }
  }
  m = m[!rmrow, ]
  m$mut[nchar(m$mut)==0] = "-" # Annotating deletions as "-"
  return(m)
}

mutations.runs = annotate_runs(mutations)

# Running dndscv

library(dndscv)
dndsout = dndscv(mutations.runs, refdb = "RefCDS_MN908947.3.peptides.rda", sm = "12r_3w", max_muts_per_gene_per_sample = Inf, max_coding_muts_per_sample = Inf, kc = "")

The code above does several things:

It identifies the sites that are different to the reference and saves them in a table of mutations.
It then collapses identical mutations in the two Omicron genomes. dNdScv treats each entry in the table of mutations as an independent mutation, so you want to make sure to input only unique mutations. If a mutation has occurred independently in two lineages, you should include it twice.
It then annotates runs of consecutive mutations into single events (again, it is important that each entry in the mutation table fed to dNdScv reflects an independent mutation).

Once you have the table of mutations, running dNdScv is straightforward, as shown in the last line of the code above. Here I used the 12r3w substitution model, which uses 12 substitution rates for all possible single nucleotide changes (equivalent to the UNREST time-irreversible substitution model).

# Saving some output tables
write.table(dndsout$sel_cv, file = "gene_selection.txt", col.names = T, row.names = F, sep = "\t", quote = F)
write.table(dndsout$annotmuts, file = "annotated_mutations.txt", col.names = T, row.names = F, sep = "\t", quote = F)
write.table(dndsout$globaldnds[c("wmis","wnon","wall"), ], file = "Genomewide_dnds.txt", col.names = T, row.names = F, sep = "\t", quote = F)

You can now look at the output of dNdScv, including:

The global (genome-wide) dN/dS ratios for missense and nonsense mutations:
Genomewide_dnds.txt
Gene-level selection analysis: This reveals strong positive selection in the spike (S gene):
gene_selection.txt
List of annotated mutations with protein-coding impact:
annotated_mutations.txt

You may be able to adapt this example to run some of your analyses. If you really do not want to use the reference genome provided and prefer to use a different reference genome, you will need to make your own CDS table and run buildref to generate a new RefCDS database. This should be easy to do for BS000685.1. You can do it manually (or automatically) using the CDS positions here. I used peptide positions for ORF1ab to have higher resolution (you can use "mat_peptide" positions from the BS000685.1 link).

I hope this helps,
Inigo

Running dNdScv on SARS-CoV-2 data #2

Description

Activity

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions