Running dNdScv on SARS-CoV-2 data

[sakrnference]

Hi
I hope you are doing well.
I can't find the full CDS table transcript for SARS-cov-2. The only one you kindly have available is for chr.name: MN908947.3.
Can you please provide me with the full trabscript table for all the chromosomes available please ?

Thank you for your time.

[im3sanger]

Hello,

Thank you. Could you clarify what you are looking for? The RefCDS object provided for SARS-CoV-2 uses MN908947.3 as the reference genome (Wuhan). The genome only has one chromosome. You can run dNdScv on SARS-CoV-2 data mapping the mutations to the MN908947.3 reference genome.

Best wishes,
Inigo

[sakrnference]

Hello,

Thank you so much for your kind response. I really appreciate your kind help.

So I want to run it on the other SARS-Cov-2 genomes references from the FASTA file I have. For example, I want to generate the RefCDS object for SARS-CoV-2 that uses BS000685.1 and other chromosome names that exists in the FASTA as the reference genomes. Was just wondering if there’s an existing CDS table transcript for that, not just for MN908947.3 ? or I am misunderstanding something here?
Also, I will really appreciate if you can please provide (if there is) or explain how to the generate the mutation list to run the dNdScv ?

I hope that was clear. Thank you so much for your time and for your help. I really appreciate it

[im3sanger]

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:

1. It identifies the sites that are different to the reference and saves them in a table of mutations.
2. 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.
3. 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:

1. The global (genome-wide) dN/dS ratios for missense and nonsense mutations:
   Genomewide_dnds.txt
2. Gene-level selection analysis: This reveals strong positive selection in the spike (S gene):
   gene_selection.txt
3. 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

[sakrnference]

Hello,

Thank you so much for your kind explanation. That was very very helpful. I really appreciate it. I have one more question I am having hard time finding the "cds.start | cds.end | length" columns values in the CDS positions link you kindly sent me. Can you please point me where can I get these values ?
Thank you so much I really appreciate your help

[im3sanger]

Hello,

I attach here the transcript table for the MN908947.3 reference (Wuhan), which shows what you need to input in those columns:
Reftable_MN908947.3_peptide_resolution-1.txt

cds.start and cds.end refer to the first and last coding base in this gene. They are required when working with CDS with multiple exones, where each exon occupies a line in the transcript table. For SARS-CoV-2 these fields are very simple. cds.start is always 1 and cds.end is the length of the coding segment (=chr.coding.end-chr.coding.start+1). cds.length is also the same length value.

For example, the S gene in the MN908947.3 reference of SARS-CoV-2 has coordinates: 21563 to 25384. cds.start is 1, cds.end and length are 25384-21563+1 = 3822. Hence the full line for the S gene in the table attached above is:

gene.id	gene.name	cds.id	chr	chr.coding.start	chr.coding.end	cds.start	cds.end	length	strand
43740568	S	GU280_gp02, spike glycoprotein	MN908947.3	21563	25384	1	3822	3822	1

Best,
Inigo

[sakrnference]

Thank you so much for your kind response. That is very helpful. So I followed your steps and I built a transcript table for the BS000685.1 and the corresponding FASTA file too but still getting the error below. I am attaching the Reftable for BS000685.1 reference and as well here is the FASTA reference for it here.
Reftable_ BS000685.1.txt

FASTA output

`DNAStringSet object of length 1:
    width seq                                               names               
[1] 29782 AGATCTGTTCTCTAAACGAACTT...AATTAATTTTAGTAGTGCTATCC BS000685.1 Severe...`

The Error I'm getting

`Error in buildref(cdsfile = cds, genomefile = fasta.data.path_test, outfile = "ouptut_refcds_test2.rda", : No chromosome names in common between the genome file and the CDS table
Traceback:`

can you please help me what is wrong here ?

thank you so much for all your help I really appreciate it so much

[sakrnference]

Hi,

I am sorry for the disturbance. I just haven't heard back from you regarding the last issue.
I would really appreciate if you kindly can help here.
Thank you so much

[im3sanger]

Hello,

Sorry for the late response, I was on holidays last week.

The error message printed by buildref provides an idea of what is causing this problem: "No chromosome names in common between the genome file and the CDS table". To build your RefCDS object, buildref needs a fasta file and a transcript table. You need to make sure that your transcript table uses the same chromosome names as the fasta file.

To solve your problem faster, please send me an email and attach your fasta file and your transcript table. If there are lessons useful for other users I can share them here.

Best wishes,
Inigo

[sakrnference]

Hello,

Thank you so much for your kind response. No worries at all I totally understand.
I double checked and both the transcript table and FASTA file has the same chromosome name.
I replied to the email you kindly sent with the files. Please let me know when you get them.

thank you so much I really appreciate your help.

[sakrnference]

Hello, Attached: transcript table and FASTA file. Please let me know how can we fix it. Thank you so much

…

On Mar 8, 2022, at 6:28 AM, Inigo Martincorena ***@***.***> wrote: Hello, Sorry for the late response, I was on holidays last week. The error message printed by buildref provides an idea of what is causing this problem: "No chromosome names in common between the genome file and the CDS table". To build your RefCDS object, buildref needs a fasta file and a transcript table. You need to make sure that your transcript table uses the same chromosome names as the fasta file. To solve your problem faster, please send me an email and attach your fasta file and your transcript table. If there are lessons useful for other users I can share them here. Best wishes, Inigo — Reply to this email directly, view it on GitHub <#2 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AU3PXMPM7QGXGXTRCKUEHUTU642U3ANCNFSM5O3JGKRQ>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you authored the thread.

-- CONFIDENTIALITY NOTICE: This E-Mail is intended only for the individual or entity to which it is addressed and may contain information that is confidential. If you have received this communication in error, do not distribute. Please delete the original message and immediately notify the sender. gene.id gene.name cds.id chr chr.coding.start chr.coding.end cds.start cds.end length strand leader-protein leader-protein leader-protein BS000685.1 212 751 1 540 540 1 nsp2 nsp2 nsp2 BS000685.1 752 2665 1 1914 1914 1 nsp3 nsp3 nsp3 BS000685.1 2666 8500 1 5835 5835 1 nsp4 nsp4 nsp4 BS000685.1 8501 10000 1 1500 1500 1 3C-like-proteinase 3C-like-proteinase 3C-like-proteinase BS000685.1 10001 10918 1 918 918 1 nsp6 nsp6 nsp6 BS000685.1 10919 11788 1 870 870 1 nsp7 nsp7 nsp7 BS000685.1 11789 12037 1 249 249 1 nsp8 nsp8 nsp8 BS000685.1 12038 12631 1 594 594 1 nsp9 nsp9 nsp9 BS000685.1 12632 12970 1 339 339 1 nsp10 nsp10 nsp10 BS000685.1 12971 13387 1 417 417 1 nsp11 nsp11 nsp11 BS000685.1 13388 13426 1 39 39 1 RNA-dependent-RNA-polymerase RNA-dependent-RNA-polymerase RNA-dependent-RNA-polymerase BS000685.1 13388 16182 1 2795 2769 1 helicase helicase helicase BS000685.1 16183 17985 1 1803 1803 1 3'-to-5'-exonuclease 3'-to-5'-exonuclease 3'-to-5'-exonuclease BS000685.1 17986 19566 1 1581 1581 1 endoRNAse endoRNAse endoRNAse BS000685.1 19567 20604 1 1038 1038 1 2'-O-ribose-methyltransferase 2'-O-ribose-methyltransferase 2'-O-ribose-methyltransferase BS000685.1 20605 21498 1 894 894 1 43740568 S GU280_gp02, spike glycoprotein BS000685.1 21509 25330 1 3822 3822 1 43740569 ORF3a GU280_gp03 BS000685.1 25339 26166 1 828 828 1 43740570 E GU280_gp04 BS000685.1 26191 26418 1 228 228 1 43740571 M GU280_gp05 BS000685.1 26469 27137 1 669 669 1 43740572 ORF6 GU280_gp06 BS000685.1 27148 27333 1 186 186 1 43740573 ORF7a GU280_gp07 BS000685.1 27340 27705 1 366 366 1 43740577 ORF8 GU280_gp09 BS000685.1 27840 28205 1 366 366 1 43740575 N GU280_gp10 BS000685.1 28220 29479 1 1260 1260 1 43740576 ORF10 GU280_gp11 BS000685.1 29504 29620 1 117 117 1

[im3sanger]

Thank you. Please send me an email with the fasta file and your transcript table to im3 at sanger.ac.uk. Link also available through my Sanger website.