text:
Plot phylogenetic tree from list of edges
I have a dataset that I wish to convert to tree or phylo format like in the ape package, in order to plot the phylogenetic tree. It is formatted like a list of edges and species names for the leaves. The dataset is quite large, but I have made an example dataset that represents it:
df <- data.frame(c(8,8,7,6,6,7,9,9,9),
                c(1,2,3,4,5,6,7,8,9),
                c(9, 9, 12, 14, 17, 19, 21, 23, 9),
                c("a","b","c","d","e", "node1", "node2", "node3", "node4"))
names(df) <- c("parent","node","branch.length", "label")

I wish to plot it as a rooted phylogenetic tree using the ape package, but I can't figure out how to convert it to for example a newick string. I have tried to look at solutions both in R and Python, but I can't figure it out. Does anyone have any idea how this can be done?
Additional information: I also have two other vectors, one contains the length of each edge, and one containing the age of the child node, equal to the sum of all edge length from that child to any leaf of the tree. It would be great if this could be included, but not neccessary for the plotting as it is only illustrative.
Thanks!
Update
I am using R version 4.1.2.
One of my attempts at translating my data to Newick code (as i understand this could include the branch lengths), using one of [these answers].
I don't understand what goes wrong, but I have not worked with Newick code before so I am not very experienced with this.
(https://biology.stackexchange.com/questions/37057/construct-phylogeny-form-edgelist):
df <- data.frame(c(8,8,7,6,6,7,9,9,9),
                c(1,2,3,4,5,6,7,8,9),
                c(9, 9, 12, 14, 17, 19, 21, 23, 9),
                c("a","b","c","d","e", "node1", "node2", "node3", "node4"))
names(df) <- c("parent","node","branch.length", "label")
tree_data <- df
str(tree_data)

tree_data$newick = paste(tree_data$label,':',signif(tree_data$branch.length,3),sep="")
nparents = 2
  while(nparents > 1){
  leaf = setdiff(tree_data$node,tree_data$parent)
  parent_leaf = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
  index_remove = integer(0)
  for(i in 1:length(parent_leaf)){
parent_children = tree_data$node[which(tree_data$parent == parent_leaf[i])]
all_children_in_leaf = all(parent_children %in% leaf)
if(all_children_in_leaf==T){
  index_remove = c(index_remove,match(parent_children,tree_data$node))
  tree_data$newick[which(tree_data$node==parent_leaf[i])] = paste("(",paste(tree_data$newick[match(parent_children,tree_data$node)],collapse=" , "),")",tree_data$newick[which(tree_data$node==parent_leaf[i])] ,sep = " ")
}
  }
  tree_data = tree_data[-index_remove,]
  parent_next = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
  nparents = length(parent_next)
  if (nparents==1) newick_code = paste("(",paste(tree_data$newick,collapse=","),")",colnames(tree_data)[parent_next],";",sep = " ")
}

tree <- ape::read.tree(text = newick_code)

plot.phylo(tree)

I don't get an error any more, but I get some unexpected output:

This is what I would expect, all nodes that are not single letters are inner nodes, not leaves:

I really appreciate any help!

There is an answer to more or less the same question in biology.stackexchange.com.

Thank you for the answer, I have taken a look at the answers but none of them seem to work for me.  (all return errors after reformatting to the data format they used) Nikunj Goel's answer of writing a code to translate to Newick code and include the branch lengths is what I have tried to do myself, but it is challenging as I am not very familiar with it.

Hi @Helle, can you provide a few lines of codes as well as the errors you are getting.   Include the version of R you are using. Those answers were published 7 yrs ago most likely you will have to adjust your code a bit to make it compatible with current standards.

@Macintosh I have updated it now :-)

the error comes from colnames(Adj)[parent_next]: Adj is not defined and should refer to the correct table.

Thanks, this helps. Now, it seems to have some trouble with the root node and "node2", they and up as leaves when I plot the tree.

yes there is indeed something off with that function. I suggest you make a comment in the post of NikunjGoel to add for advice! You can update the question with the Newick output and the expected newick output

Unfortunately I don't have enough reputation comment or message anyone!

I have added a comment to that answer.

The following example uses python to generate the newick string. You can than import that in R for further analysis.
Please keep in mind that you can not use self reference for node4 (9,9,9), so I changed this to 'Root'.
from ete3 import Tree
import pandas as pd

df = pd.DataFrame({'parent': [8, 8, 7, 6, 6, 7, 9, 9, 'Root'],
                   'node': [1, 2, 3, 4, 5, 6, 7, 8, 9],
                   'branch.length': [9, 9, 12, 14, 17, 19, 21, 23, 9],
                   'label': ['a', 'b', 'c', 'd', 'e', 'node1', 'node2', 'node3', 'node4']})

parent_child_list = df[['parent', 'node', 'branch.length']].apply(tuple, axis=1).tolist()
tree = Tree.from_parent_child_table(parent_child_list)
newick_str = tree.write(format=1)

newick_str:
'(((3:12,(4:14,5:17)6:19)7:21,(1:9,2:9)8:23)9:9);'

Thank you! Exactly what I was looking for.


text:
Which human reference genome version does VARAdb use?
I am interested in annotating a list of variants using VARAdb. My coordinates use the GRCh38 (hg19) reference genome. However, I do not know which coordinate system VARAdb uses.
A section of the VARAdb publication reads:
"GWAS have provided a large amount of data
associating genetic variants with common phenotypes. We
collected risk SNPs fromfive sources, including theNHGRI
GWAS Catalog (29), GWASdb v2.0 (30), GAD (31), Johnson
and O’Donnell (32), and GRASP v2.0 (33). Then, we
filtered risk SNPs ofwhich ‘Variant ID’which are not ‘rsID’
and transformed locations of risk SNPs from the GRCh38
assembly version into the hg19 version in the GWAS Catalog.
Finally, we obtained 1 515 001 risk SNPs associated
with diseases, traits and phenotypes."
which makes me think that they may be using the hg19 genome, but I am not sure.
So, my question is:
Which human reference genome version does VARAdb use?
VARAdb database:
http://www.licpathway.net/VARAdb/search.php
VARAdb publication:
https://pubmed.ncbi.nlm.nih.gov/33095866/

Good you've got the answer. FYI @CristianRiccio you'll find you can now upvote because you've >20 rep  (hint, hint :-) ).

They are using hg19.
I found this by simply going to "Browse" then copying taking the first rsID in the list (rs1891805) and comparing the chromosomal position they give in VARAdb to what I see in dbSNP. VARAdb gives it as position 154904840 on chromosome 1 and dbSNP has:
Sequence name       Change
GRCh37.p13 chr 1    NC_000001.10:g.154904840T>C
GRCh38.p14 chr 1    NC_000001.11:g.154932364T>C
PMVK RefSeqGene     NG_053028.1:g.15287A>G

GRCh37 is hg19, and GRCh38 is hg38. So, as you can see above, the position 154904840 is from hg19 and the hg38 position for this rsID is 154932364.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
What tool(s) can I use to use a .faa file to feed into JSMol?
JSMol, the molecular modelling software, requires PDB file format to work and render the graphics. Currently, I am writing an annotation software - backend in Perl - that takes in a .faa file (fasta amino acids). Is there any way or tool out there that can easily convert to a PDB file? Many of the tools out there do not seem to give what I need.
My workflow so far:
input.fna -> Bakta -> bakta.faa -> covert tool -> proteins.pdb -> JSMol
Essentially, a tool exists called pdb2fasta that takes a pdb file and converts to a fasta file but can’t seem to find a proper tool that does it the other way around.

Or do you just want an unfolded sequence of amino acids (as atoms and bonds)? In any case, you will need something that generates those coordinates, even if just as a linear 'model' of a chain

Would you happen to have any ideas what I could use to do that? @gilleain

@HarishPrabhakar That is an exceptionally difficult problem for a variety of reasons. there is not an obvious way to translate linear AA sequence into three-dimensional structure. You may have heard of alphafold, which uses neural nets to more or less solve the problem. But you have to understand that this is a stochastic and error-prone prediction!

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to retrieve sequences from a txt file by protein names from other txt file?
My problem is that I want to retrieve the sequences line + secondary structure line that match the protein names from the header.txt file:
The text file(header.txt) contains the protein names that look like this:
1XXIJ
1P0EA
1H9HI

And the other file(pdb.txt) structured like this:
Line 1:the protein name
Line 2: protein length
Line 3: amino acid sequence
Line 4: secondary structure
Line 5: empty
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

So my expected output is:
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

Do you guys have any suggestions for me by awk or any(Linux)code? Please let me know if you require additional information. Cheers!

If your file structure is consistent, one potential solution is to use GNU grep:
grep -A4 --no-group-separator -Fwf header.txt pdb.txt

This works with your example data and should be relatively fast on your 'real' data, due to -Fw ('fixed' i.e. not regex, and only match whole words) but it requires GNU grep for the --no-group-separator option (i.e. not BSD grep).

The only issue here is that it tests $1 in a for each line of input (in pdb.txt). It might not do the right thing if the list of headers contains a protein name called 'M', for example. It's exactly the same as: FNR==NR { a[$1]; next } $1 in a { x = FNR+4 } FNR <= x.

I'm assuming that your second file is tab-delimited and has no commas. The easiest way to do this would be to transpose each record in your second file to a single line and then grep, then transform the result again:
grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"

I just tested and it works:
➜ cat pdb.txt
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

➜ cat header.txt
1XXIJ
1P0EA
1H9HI

➜ grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S
```

One way using awk:
awk 'FNR==NR { a[$1]; next } FNR%5==1 { f = ($1 in a) ? 1 : 0 }f' header.txt pdb.txt

This adds each of the protein names in header.txt to an array. By adding the first field ($1), rather than the complete record ($0) to the array, we can ensure that we strip off any lagging whitespace which can sometimes cause problems depending on how the list of protein names was created. Next, we can check whether the record number (i.e. line number) of the second input file is a multiple of 5 with remainder 1. This condition returns true for the first line of each record in pdb.txt. We can then use a ternary operator to set a flag variable to one if the first field of the current line (i.e., the protein name) is present in the a array, or to zero to indicate that the current record is not one of the records in the array. The f by itself is just shorthand for f != 0 { print }, which prints the current line when the flag is non-zero.

Thanx guys. All of the codes above work well for me. I appreciated it.


text:
Why does this nextflow script finish after running one sample?
I have 36 samples in total in the bam_files folder, with the name like this 20230306_CH_EP_C01.md.bam
I expect the code to output all the 36 samples, one by one. But the run stopped after only running one sample.
I have searched online, but did not find a solution.
I am using nextflow version 2.
bam_chl = Channel
     .fromPath("bam_files/*.md.bam")
     .map{ item -> tuple(item.baseName.substring(0, 18), item) }

ref_chl = Channel
      .fromPath("sequences/hg38.fa")
ref_index_chl = Channel
      .fromPath("sequences/hg38.fa.fai")
amplicon_chl = Channel
      .fromPath("sequences/tp53_amplicon.hg38.interval_list")
target_chl = Channel
      .fromPath("sequences/TP53.hg38.interval_list")

process get_pcr_metric {

  publishDir path: 'pcr_metric', pattern: '*.csv'

  input:
  tuple val(pair_id), path(bam)
  path(ref)
  path(ref_index)
  path(amplicon)
  path(target)

  output:
  path("*.pcr_metrics.csv")

  """
  java -jar /picard.jar CollectTargetedPcrMetrics \
       I=$bam \
   O=${pair_id}_pcr_metrics.txt \
       R=$ref \
   AMPLICON_INTERVALS=$amplicon \
       TARGET_INTERVALS=$target

  sed -n -e 7p -e 8p ${pair_id}_pcr_metrics.txt |  sed 's/\t/,/g'   > ${pair_id}.pcr_metrics.csv

  """

}

workflow {

get_pcr_metric(bam_chl, ref_chl,ref_index_chl, amplicon_chl, target_chl)
}

The issue here is that the fromPath factory method produces a queue channel, but you don't want all of your input channels to be queue channels. Most of the time, what you want is one queue channel and one or more value channels when you require multiple input channels.
You could use the value factory method to create the value channels, but a simpler way is to just pass in some file objects. The latter works because a value channel is implicitly created by a process when it is invoked with a simple value1. For example:
params.ref_fasta = "./sequences/hg38.fa"
params.amplicon_intervals = "./sequences/tp53_amplicon_intervals.hg38.interval_list"
params.target_intervals = "./sequences/TP53.hg38.interval_list"

params.bam_files = "./bam_files/*.md.bam"
params.outdir = './results'

process get_pcr_metric {

    publishDir path: "${params.outdir}/pcr_metric", mode: 'copy'

    input:
    tuple val(pair_id), path(bam_file)
    path ref_fasta
    path amplicon_intervals
    path target_intervals

    output:
    path "${pair_id}_pcr_metrics.txt"

    """
    picard CollectTargetedPcrMetrics \\
         I="${bam_file}" \\
         O="${pair_id}_pcr_metrics.txt" \\
         R="${ref_fasta}" \\
         AMPLICON_INTERVALS="${amplicon_intervals}" \\
         TARGET_INTERVALS="${target_intervals}"
    """
}

workflow {

    ref_fasta = file( params.ref_fasta )
    amplicon_intervals = file( params.amplicon_intervals )
    target_intervals = file( params.target_intervals )

    Channel
        .fromPath( params.bam_files )
        .map { tuple( it.getBaseName(2), it ) }
        .set { bam_files }

    get_pcr_metric( bam_files, ref_fasta, amplicon_intervals, target_intervals )
}

The code can be made to work on all samples by changing the interface of the process and using the combine operator:
ref_chl = Channel
      .fromPath("sequences/hg38.fa")
ref_index_chl = Channel
      .fromPath("sequences/hg38.fa.fai")
amplicon_chl = Channel
      .fromPath("sequences/tp53_amplicon.hg38.interval_list")
target_chl = Channel
      .fromPath("sequences/TP53.hg38.interval_list")

bam_chl = Channel
     .fromPath("bam_files/*.md.bam")
     .map{ item -> tuple(item.baseName.substring(0, 18), item) }
     .combine(ref_chl)
     .combine(ref_index_chl)
     .combine(amplicon_chl)
     .combine(target_chl)

process get_pcr_metric {

  publishDir path: 'pcr_metric', pattern: '*.csv'

  input:
  tuple val(pair_id), path(bam), path(ref), path(ref_index), path(amplicon), path(target)

  output:
  path("*.pcr_metrics.csv")

  """
  java -jar /picard.jar CollectTargetedPcrMetrics \
       I=$bam \
   O=${pair_id}_pcr_metrics.txt \
       R=$ref \
   AMPLICON_INTERVALS=$amplicon \
       TARGET_INTERVALS=$target

  sed -n -e 7p -e 8p ${pair_id}_pcr_metrics.txt |  sed 's/\t/,/g'   > ${pair_id}.pcr_metrics.csv
  """
}

workflow {
  get_pcr_metric(bam_chl)
}

The combine operator as the name implies generates all combinations. The different file objects are then in one tuple.
Since only one Channel is passed to the process, it's now also possible to rewrite the workflow code as follows:
workflow {

  refs = [
    "sequences/hg38.fa",
    "sequences/hg38.fa.fai",
    "sequences/tp53_amplicon.hg38.interval_list",
    "sequences/TP53.hg38.interval_list"
  ]

  ref_chl = Channel.from(refs)
   | map{ pointer -> file(pointer) }
   | collect

  Channel.fromPath("bam_files/*.md.bam")
   | map{ item -> tuple(item.baseName.substring(0, 18), item) }
   | combine( ref_chl )
   | get_pcr_metric

}

You do not need the combine operator here: In general, multiple input channels should be used to process combinations of different inputs, using the each qualifier or value channels.


text:
Gene symbol list for all protein coding genes in mice
How can I get a csv or a list of the gene symbol names for all the protein coding genes in mice?  I have RNA sequencing data and I'm not interested in the non-coding stuff. I'm worried it could mess up my results.  Can GEO or blast provide this?
Edit: I am not interested in EntrezID or Ensemble, I don't not want them I just am more interested in gene name like Fos or Penk. Etc.  I work in RNA seq I don't need to map any fastq files to a genome that has been done already I just need to filter out genes that are not protien coding.  If you are not familiar with what RNA seq data looks like it is simply rows of genes and columns of samples each entry giving the number of reads for that gene in that sample, single cell RNA seq is the same.  Have a look at GEO if you want more examples: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE145970

So my answer should give you exactly what you need, right? If you don't want protein coding genes that do not have a "name" assigned, you can skip the "gene stable id" option.

If you want all known mouse protein coding genes, you can get the list from Ensembl's BioMart.
First, select the mouse genome:

Then, in "Filters" limit to protein coding genes only:

And next, choose "Gene name" and, to have a more standardized and stable name and to ensure that you get allr esults (not all genes have a "name", but they will all have IDs), include the "Gene stable ID":

Finally, click on "Results":

Select "File" and "CSV" and then click on "Go" to download the list:

The end result is a file like this:
Gene stable ID,Gene name
ENSMUSG00000064341,mt-Nd1
ENSMUSG00000064345,mt-Nd2
ENSMUSG00000064351,mt-Co1
ENSMUSG00000064354,mt-Co2
ENSMUSG00000064356,mt-Atp8
ENSMUSG00000064357,mt-Atp6
ENSMUSG00000064358,mt-Co3
ENSMUSG00000064360,mt-Nd3
ENSMUSG00000065947,mt-Nd4l
ENSMUSG00000064363,mt-Nd4
ENSMUSG00000064367,mt-Nd5
ENSMUSG00000064368,mt-Nd6
ENSMUSG00000064370,mt-Cytb
ENSMUSG00000095019,
ENSMUSG00000094915,
ENSMUSG00000095041,
ENSMUSG00000095666,
ENSMUSG00000094741,Pramel39

With 21849 gene IDs, of which 21772 have a "name" as well.

Thanks, @RamRS, I blame my years living in France! :)


text:
What is the best approach to do a PCA analysis to distinguish between two groups?
I have two samples A and B and I do have an RNAseq for both of them. Inside each file (for A and B) I do have a set of FPKMs that give me a fold change for each gene (sample vs control). It looks somewhat like that:
gen_name grp01 grp02 grp03 grp04 fold_change
gene1    10    10    5     5     2
gene2    11    11    11    11    1
gene3    5     5     10    10    0.5
gene4    ...

for sample A
and, say
gen_name grp05 grp06 grp07 grp08 fold_change
gene1    100   100   5     5     20
gene2    110   110   110   110   1
gene3    5     5     100   100   0.05
gene4    ...

for sample B
I need to show with PCA that, foe example, gene1 is sample B is much more overexpressed than in sample A, same for gene3, but it looks like gene2 is not very interesting.
What data must be put in PCA analysis?

Using a package is the best approach. Scanpy is popular right now. I do this stuff via first principles, so my advice isn't cool here. Anyway, I'd combine all data into one dataframe, minus log-fold change. I assume log fold change is the opposing sample? The data must have a column to identify A genes from B genes, or denoted within the gene name. PCA does not give gene by gene comparisons, it will be cluster by cluster. Interpreting each clusters is done independently (by mapping the log-fold-change), i.e. what is clustering together and that cluster represents a cohesive expression group.

The alternative is DEG analysis, gene by gene - which you've already done.

Thank you. I want to show that A is different from B based on an expression of gene1 and gene3 (in this test experiment). I want to do this in R but scanpy is a python package. Is there any way to use prcomp() ?


text:
Identifying somatic mutations in cell lines
I would like to identify the somatic mutations present in a cell line and characterise the genes that are potentially affected by those mutations. For example, are there oncogenes mutated in a subpopulation in the cell line?
I have currently access to single cell WGS data and sequencing data of the bulk sample derived from the same clone. The data was generate after MDA amplification of the single cell's genome so I am currently evaluating SCCaller and Prosolo for variant calling. If I correctly understand single cell variant calling, mutations called in each cell are called relatively to the bulk background sample. So each cell separately doesn't provide much information on the cell line itself.
My questions are:

Do i need single cell WGS to quantify the overall 'genomic health' of a cell line or can I just use deep sequencing followed by somatic variant calling against a reference genome.
What added value can single cell WGS provide here and how can I extract that information?

Links to papers, packages or statistical models more than welcome.
Thanks!

I am not actively following this field, so take my answer with a grain of salt.
With bulk data, you can't easily tell how subclones differ from each other or characterize the mutation rates in subclones. It is also challenging to get accurate somatic mutation calls especially when you don't have the matched normal.
With single-cell data, you usually get more accurate point mutation calls in 20-80% of the genome, depending on the amplification method. You can estimate the mutation rate and call large CNVs more accurately in each cell. Sometimes you can even reconstruct the evolution of subclones. I would say single-cell data is better than bulk data all around. The only problem is that the single-cell approach is more costly as you need to sequence many cells to see the pattern.
On a side note, MDA is not that accurate for calling point mutations. This may not be a problem if the mutation rate is much higher than the error rate.


text:
Converting a Dante Lab VCF file
Is there an easy way to grab the rsids (hg38) of a VCF file? I know one of the tabs can contain RSIDs, but this file doesn't contain any.
I have a software I created for mapping 23AndMe DTC testing to SNPedia, and I'm trying to integrate more DTC testing. (https://www.github.com/mentatpsi/osgenome)
The user found this https://gist.github.com/KlausTrainer/26c2996cf32677e4c107bd8aaee67794 which is a script to use variant summary data to map to the RSIDs, but I'm hesitant to perform a SQL dump to find a table that contains
chromosome, postion, and rsID columns.
As a bonus question, 23AndMe reports all positive orientations. How would I go about converting strands to match this orientation? Does the VCF contain information that would allow this action autonomously?

Can you show us an example of your file? VCFs do not need to contain rsIDs and often don't, as you see. Also, if you are referring to the raw, unfiltered data that you can download from 23&me, you need to be very, very careful with them because, at least when I looked at them last, they don't contain any quality information and will contain many, many false positives. Which is why they are not included in the official results you get on their webpage.

I do have a disclaimer on my application that Raw Data coming from DTC vendors have a false positive rate approaching 40% per a research article found in a Bioinformatics magazine.

If my understanding of the data is correct... this is the format: https://www.beholdgenealogy.com/blog/?p=2879

Re the false positives - with this VCF data, there appears to be several fields of quality data. How useful I don't know, being a beginner to this subject. I'll post a  the header and a few lines of the VCF later in a reply. But if you have, for instance, the number of reads (is that depth? Still getting the terminology) that match this particular variant, I assume that would allow you to provide a confidence level.

If your file doesn't contain the rsIDs, you can use a tool like Ensembl's VEP to annotate the file and get rsIDs as well. Of course, not all variants have an rsID so don't expect to find one for everything reported in your file.
More importantly, if these are the raw, unfiltered data from commercial direct-to-consumer sequencing products like 23&me, you need to be aware that they are incredibly untrustworthy. This is why 23&me don't include them in their official report. They only show you the quality results that have been filtered. When working with the raw data, you are very likely to have completely bogus results and will need to have a full VCF file and filter on various quality criteria.
I have looked at the raw 23&me data from 4 individuals, and all of them had extremely rare pathogenic variants reported that would basically have killed them long before they learned to walk. They were clearly false positives.
I stress that I am not suggesting these companies are producing crap: none of these variants was reported in the official report of 23&me, they were filtered out. The bad calls are only present when you work on the raw, unfiltered data but if that is what you are working on, you  need to be extremely careful.

Is there a downloadable file which 'maps' chromosome, position to rsid? According to this https://bioinformatics.stackexchange.com/a/14193/17155 one can download a dbSNP database file, with caveats - does this file allow you to map chr:pos to rsid? The caveats mentioned in that answer are that dbSNP build numbers are important (range mentioned is 152, 153, 154) but that reference genome is not. What build is suitable for use with SNPedia? Would a suitable file for Dante, using GRCh37, be say https://ftp.ncbi.nih.gov/snp/pre_build152/organisms/human_9606_b150_GRCh37p13/VCF/00-All.vcf.gz

Ensembl's VEP, that you mentioned, seems to return vast amounts of data, tens of rows per chr:pos, with reams of redundant data, and a limited query size, and very long waits even for tiny queries. Unless I am misunderstanding and the 'mapping' is a lot more complicated than chr:pos to rsid, it will be easier to just do it locally on my PC after a longish single download of a few GB aided by a bit of fast C code and an SSD.

@Pete are you really only giving chr and position, with no actual variant? You said you had a VCF file, do you not?

Sorry, you're absolutely right, I have a long line of data associated with each chr:pos - format is chr, pos, <period (in place of rsid)>, reference base, my actual variant base(s), and then there's several quality fields and many more fields. The fields are described in the file header, there must be about 30 more fields iirc (I'm not in front of the file right now). Many of these fields I don't yet understand. One thing I'm still unsure of is if you have 2 different possible alternative bases at a particular position, do they get different rsids or is it the same rsid with different 'types'

@Pete That... isn't a VCF and isn't even similar to the image of a file you showed. That's why I asked you for an example of your file. Instead, you gave me a photograph (images of text are useless, we cannot run text parsing tools on them) and now you tell me it was a photograph that isn't even related to your data? Please don't do this. That invalidates the entire question. So please go and edit your question and include some actual lines from your file and the output you expect from them. As for rsIDs, sometimes one rsID corresponds to a specific variant, other times it can be many.

No the above linked blog post is not me, but I'm the guy who contacted Shay originally about his app, and he was looking into it for me. Anyway I'll post a section of my own file in a few minutes, need to get to a pc.


text:
How to tackle an Index Error in my code for finding ORFs
I have created a python code that prints out all lengths of ORFs in a DNA sequence (start codon - "ATG", stop codons - "TAG", "TAA" and "TGA". It seems to give correct answers, when compared with specialized ORF-finding programs. The thing is that besides correct answers it also gives an Index Error in the end:
Traceback (most recent call last):
  File "scripts.py", line 33, in <module>
    if list_start[i]>list_stop[n-1]:
IndexError: list index out of range

Now, I have created a dictionary from a FASTA file where key is a sequence name and value is the sequence itself and want to find ORFs for all sequences in that dictionary using a for loop and this code. However, the loop runs only once and I think that the Index Error to be the cause.
The code is the following:
#!/usr/bin/env python3
line="TCGGGCGAAGGCGGCAGCAAGTCGTCCACGCGCAGCGCGGCACCGCGGGCCTCTGCCGTGCGCTGCTTGGCCATGGCCTCCAGCGCACCGATCGGATCAAAGCCGCTGAAGCCTTCGCGCATCAGGCGGCCATAGTTGGCGCCAGTGACCGTACCAACCGCCTTGATGCGGCGCTCGGTCATCGCTGCATTGATCGAGTAGCCACCGCCGCCGCAAATGCCCAGCACGCCAATGCGTTCTTCATCCACATAGGGGAGCGTTACGAGGTAGTCGCAGACCACGCGGAAATCCTCGACGCGCAGTGTCGGGTCTTCGGTAAAACGTGGTTCGCCGCCGCTGGCACCCTGGAAGCTGGCGTCGAAGGCGATGACGACGAAACCTTCCTTGGCCAGCGCCTCGCCATACACGTTCCCCGATGTTTGCTCCTTGCAGCTGCCGATCGGATGCGCGCTGATGATGGCGGGATATTTCTTGCCTTCGTCGAAGTTCGGCGGGAAGTGGATGTCGGCTGCGATATCCCAATACACATTCTTGATCTTGACGCTTTTCATGACAGCTCCGTTCAGGGGGAGGGGGTAAGTTCGCCAGGCCGAATCGTTGGTAGCCAAGCGGCAACGACTCGAATATAGAGAGCCGATTGGAATTCCGTAAGATCGCAATCTGGACTACAGTGGTATCTTCAAATTGACAATGGCACCTACATGGATCCCTCACTGCTTCCGTCTCTCGCGTGGTTCGCCCACGTCGCACATCATCGTAGCTTCACGAAAGCGGCTGCGGAAATGGGCGTTTCTCGAGCAAACCTGTCGCAGAACGTGAAGGCGCTCGAACGCCGGTTGAACGTCAAGCTGCTGTATCGAACGACTCGCGACATGTCGCTGACCGAGGAGGGGCAGCGGCTCTACGAGGTGTGGTATCCCGCGCTGGTCGCGGTCGAGCGGACGGTCGACGCGCTGCACGAGGAGCGCGACGAGCCGTCGGGGTTGATTC"
r_frame=input('Enter a reading frame')
line=line[int(r_frame)-1:]
# making list of codons (3-nucleotide groups)
x=3
codon=[line[y-x:y] for y in range(x, len(line)+x,x)]
#defininig function  that will find all the indexes of codon in list of codons
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(codon):
        if value == item_to_find:
            indices.append(idx)
    return indices
# makirng list with all the indexes of "ATG" (start codon)
list_start=find_indices(codon, "ATG")
#making list with all the indexes of "TGA", "TAA" or "TAG" (stop codons) and sorting it from smaller to bigger
list_stop=find_indices(codon, "TGA")+find_indices(codon, "TAA")+find_indices(codon, "TAG")
list_stop.sort()
i=0 # variable for indexes in list_start
n=0 #variable for indexes in list_stop
j=0 # variable that determines numer of times searching for ORF will be repeated
#first round without loop so that i and n will be >0 and we can make  if list_start[i]>list_stop[n-1] without errors
if list_start[i]<list_stop[n]:
    print(3*(list_stop[n]-list_start[i])+3)
    j=j+1
    n=n+1
    i=i+1
else:
    n=n+1
#loop that will search for suitable start codon-stop codon pairs and print the length of ORF
while j<codon.count('ATG'):
    if list_start[i]>list_stop[n-1]:
        if list_start[i]<list_stop[n]:
            print(3*(list_stop[n]-list_start[i])+3)
#multiply by 3 'cause one item in list stands for 3 nucleotides and add 3 so the whole length of stop codon will be counted
            j=j+1
            n=n+1
            i=i+1
        else:
            n=n+1
    else: i=i+1

For sequence given in the code the ORFs lengths should be as following:
Reading frame 1:63; 36; 36; 213; 129
Reading frame 2: 39
Reading frame 3: 186
This code might seem weird as I'm new to python and bioinformatics. I will be very thankful if anyone gives me an advice how can I get rid of the Index Error here.

Firstly, there is a bug in the code,
Original,
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(codon):
        if value == item_to_find:
            indices.append(idx)
    return indices

Corrected,
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(list_to_check):
        if value == item_to_find:
            indices.append(idx)
    return indices

You renamed the variable codon for the function, but then retained its name as codon within the function. It should be list_to_check. Also it is convenient.
If you are certain your code is fine, the easiest way to silence the error is via a try and except statement. The critical variable is line and I'm not sure where that is being declared. Without that its difficult to unpick the underlying error, which is an out of range error.
def find_indices(list_to_check, item_to_find):
        indices = []
        for idx, value in enumerate(list_to_check):
            if value == item_to_find:
                indices.append(idx)
        return indices

if __name__ == '__main__':
    #line= sample DNA nucleotide sequence

    try:
        r_frame= input('Enter a reading frame')
        line = line[int(r_frame)-1:]
        line=line[int(r_frame)-1:]
        x=3
        codon=[line[y-x:y] for y in range(x, len(line)+x,x)]
        list_start=find_indices(codon, "ATG")
        list_stop=find_indices(codon, "TGA")+find_indices(codon, "TAA")+find_indices(codon, "TAG")
        list_stop.sort()
        i=0 # variable for indexes in list_start
        n=0 #variable for indexes in list_stop
        j=0 # variable that determines numer of times searching for ORF will be repeated
        if list_start[i]<list_stop[n]:
            print(3*(list_stop[n]-list_start[i])+3)
            j=j+1
            n=n+1
            i=i+1
        else:
            n=n+1
        while j<codon.count('ATG'):
            if list_start[i]>list_stop[n-1]:
                if list_start[i]<list_stop[n]:
                    print(3*(list_stop[n]-list_start[i])+3)
                    j=j+1
                    n=n+1
                    i=i+1
                else:
                    n=n+1
            else:
                i = i + 1
    except AssertionError as error:
        print (error)
    except:
        pass

If you want to see the error you simply deleted the late two lines of code.
However* to resolve the second bug is quite complicated, particularly without a line declaration. We would simply use is Biopython's .translate() method, we wouldn't do this from first principles. In a past question thats been done.


text:
which secondary structure assignment algorithm is used by PDB?
older releases of PDB used to specify DSSP as the 2ary structure assignment algorithm used in PDB, right now no algorithm is explicitly stated. Where can I find a reference for it in the database?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Does the PDB assign the secondary structure tags, or just validate the one uploaded in the repository ??

DSSP 4.0 ?? https://github.com/PDB-REDO/dssp is used in PDB-redo

cant find any reference tough


text:
Replace a string in a column of a text file based on matching string listed in another file using linux
I have multiple sets of files, so I grepped the list from a column in a seperate file such as file1.txt. I want to change all other files by using this file1.txt based on a matching string.
for example I am pasting three files here
File 1:
chr - seq1 NZ_JAHWGH010000010.1 0 60562 green_a4
chr - seq3 NZ_JAHWGH010000012.1 0 466573 green_a4
chr - seq5 NZ_JAHWGH010000013.1 0 125526 green_a4
chr - seq39 NZ_CAJGBF010000006.1 0 211154 green_a4
chr - seq40 NZ_CAJGBF010000007.1 0 199119 green_a4
chr - seq41 NZ_CAJGBF010000008.1 0 190934 green_a4
chr - seq42 NZ_CAJGBF010000009.1 0 178251 green_a4

I extracted the required information using grep; the command is given below:
grep "NZ_JAHWGH" file.txt |awk '{print $3}' >file1.txt
grep "NZ_CAJGBF" file.txt |awk '{print $3}' >file2.txt

Now I wanted to use this file1.txt and file2.txt to change another files. Output will be like if list in file1.txt matches then all matched lists should be substituted as seq1 similarly using the second file2.txt I want to substitute all matching patterns as seq2.
For example the file is given here to be changed by using information from the previous file:
seq3 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq5 167934 169650 fill_color=green_a0
seq42 200621 202367 fill_color=red_a0
seq41 160164 161903 fill_color=green_a0
seq40 31356 33104 fill_color=green_a0

And final output after the change should be as follows:
seq1 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq1 167934 169650 fill_color=green_a0
seq2 200621 202367 fill_color=red_a0
seq2 160164 161903 fill_color=green_a0
seq2 31356 33104 fill_color=green_a0

I tried using this command:
for i in f2.txt; do grep -wf f1.txt $i|awk '{$1="seq1";print}'; grep -wf other.txt $i |awk '{$1="seq2";print}'; done

2nd attempt:
grep -wf f1.txt f2.txt|sed -e 's/seq[0-9]*/seq1/g'

Is this correct? If it is correct is there any more efficient way to do it? also it is only printing to the console. I want the output and also changes in the main file. Kindly improve my scripts.
I appreciate yout time and help.

So if I understand this correctly you have a file with multiple seqs and depending on a column value for example ""NZ_JAHWGH" you want to rename the seqx term with seq1, seq2, seq3 etc, depending on the column value. In the end all seqx values associated with NZ_y will be seq1, NZ_z will be seq_2 etc.. Is that right? It might be helpful to specify the context somewhat

ya its like that, First I want to rename first file, and later same changes with other files also, but main information is only provided in first file such as information of sequence, so similar changes should also needed in other file since every file will be used for analysis in one go. Thanks for ypur comment

Here's one way using awk. This avoids the creation of the intermediary files by matching the alphabetic characters in the fourth column. If the prefix has not been seen, it is added to an array along with a counter, which is only incremented when the prefix has not yet been seen. We then store the seq ids and the replacement strings in another array. Once we've finished reading the first file, we can then query the first column of the second file in the arguments list. If the seq id in the first column exists in the array of seq ids, we can print the replacement string and remaining columns.
With the following in a file called script.awk:
BEGIN {

   counter = 1
}

FNR==NR {
    if( match($4, /^NZ_[A-Z]+/) ) {
        prefix = substr($4, RSTART, RLENGTH)

        if( !(prefix in seen) ) {
            seen[prefix] = counter++
        }

        seq_ids[$3] = "seq" seen[prefix]
    }

    next
}

$1 in seq_ids {

    print seq_ids[$1], $2, $3, $4
}

Run using:
awk -f script.awk file1 file2

Results:
seq1 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq1 167934 169650 fill_color=green_a0
seq2 200621 202367 fill_color=red_a0
seq2 160164 161903 fill_color=green_a0
seq2 31356 33104 fill_color=green_a0

thank you very much,


text:
Project PLINK eigenvector from one group on top of another group
I have 1 .bcf file called with individuals from two groups, one which is my samples, another which is a reference panel. Due to my samples having a limited amount data, when I run PLINK to generate .eigenvec, none of the samples have any overlap. I believe this is due to the data being skewed. Therefore, I wish to generate a PCA from my reference panel and my samples separately, and then project my samples on top of the reference panel.
I tried from the PLINK file from the original .bcf file to extract the two groups, then identify which genotypes are equally well represented using parameters such as --geno 0.1 to generate SNP's overlapping between the groups. Then I tried the following after reading a lot of different suggestions online, which generates a nice plot with the points projected on top of each other.
group1_eig <- fread("Samples.eigenvec")
group1_eig <- group1_eig[,3:ncol(group1_eig)]

group1_freq <- fread("Samples.frq")
group1_eval <- fread("Samples.eigenval")

# read in the files
group2_eig <- fread("RefPanel.eigenvec")
group2_eig <- group2_eig[,3:ncol(group2_eig)]

group2_freq <- fread("RefPanel.frq")
group2_eval <- fread("RefPanel.eigenval")

freq <- rbind(group1_freq, group2_freq)
eval <- rbind(group1_eval, group2_eval)
eig <- rbind(group1_eig, group2_eig)
pca <- prcomp(eig[,3:ncol(eig)], scale = TRUE)
group1_pca <- predict(pca, newdata = group1_eig[,3:ncol(group1_eig)],
                      type = "scores")

Is this approach acceptable? What other approaches work? I have read you can also normalise one group using the square root of the eigenvalues for the other group. Or potentially use the MAF (.frq).


text:
How can I programmatically retrieve all alternative gene names from Wormbase?
I'm trying to get all C. elegans genes from Wormbase programmatically, as well as all of their alternative gene names. To do this, I'm using intermine.
I was able to obtain all genes, but not ALL alternative names. The code below retrieves the genes' primary identifiers, and their "secondary identifiers."
from intermine.webservice import Service
service = Service("http://intermine.wormbase.org/tools/wormmine/service")

query=service.new_query("Gene")
# 6239 is the taxonId for C. elegans
query.add_constraint("organism.taxonId","=","6239")
query.select("*")

all_wormbase_genes_symbol = list()
all_wormbase_genes_sID = list()

for row in query.rows(start=0):
    # get primary symbol/ID
    all_wormbase_genes_symbol.append(row["symbol"])
    # get secondary ID
    all_wormbase_genes_sID.append(row["secondaryIdentifier"])

I tried to look for other keys in row that might contain this information by first retrieving the list of all possible keys:
k = set()
for row in query.rows(start=0):
    k.update(row.keys())

but none have what I want:
{'Gene.automatedDescription',
 'Gene.biotype',
 'Gene.briefDescription',
 'Gene.description',
 'Gene.geneticPosition',
 'Gene.id',
 'Gene.lastUpdated',
 'Gene.length',
 'Gene.name',
 'Gene.operon',
 'Gene.primaryIdentifier',
 'Gene.score',
 'Gene.scoreType',
 'Gene.secondaryIdentifier',
 'Gene.symbol'}

I don't know via intermine, but you can download a file with this information from the FTP at /pub/wormbase/releases/WS287/species/c_elegans/PRJNA13758/annotation/c_elegans.PRJNA13758.WS287.geneOtherIDs.txt.gz (changing to the appropriate WSxxx version)


text:
How to generate rsid values locally from VCF file
I recently used Dante to obtain a WGS and have downloaded several files produced during this process. There is a file that ends '.vcf' which I did think was a VCF file but I'm not sure now after a previous question. Part of that file is included below. Its about 1.1GB in size.
It appears to me to follow VCF format at least to a degree but it is missing the rsid's. The field where an rsid would usually be in a VCF file is replaced with a period (.).
I would like to backfill (annotate?) this file with known rsid's (where defined). There are nearly 4 million of them and online tools don't seem ideal, I would prefer to be able to do it locally, if that is possible. I am not looking for software to do this, but a data file which specifies which rsid is assigned to a particular chr:pos:allele. That may or may not be possible - perhaps I misunderstand and its far more complex than that. But from what I can tell, for a particular reference genome, and a particular build of dbsnp, there is a simple relationship as described.
I want to do this, for instance, to be able to use various tools, like SNPedia. Dante appears to use GRCh37, vs SNPedia using mainly GRCh38, which I can see is an issue, and will probably require me to buy a GRCh38 conversion from Dante (or do it myself). But please ignore that particular issue right now.
I don't quite understand what the dependencies are. I've seen comments that suggest that reference genome is not important to generating rsid's from VCF files and that only dbsnp build number needs to match the environment that the VCF file was generated in, but that does not seem to make much sense to me. From what I can see, if you are converting chr:pos information to rsid's, you need to be using the right reference genome. Does dbsnp build really matter? Should you not use the latest build, if it is continually being added to - surely rsid's are not reused? And then a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.
So in terms of question - is there a file, downloadable from somewhere, that contains mappings from chr:pos to rsid, and if so, how does one select the correct file for a particular dataset (IE generated with a particular reference genome, with a particular dbsnp, to be used with references that have used a particular reference/dbsnp build like GRCh38 / dbsnp build 151)?
This file (header first), downloaded from Dante, is named .filtered.snp.vcf
    ##fileformat=VCFv4.2
    ##FILTER=<ID=PASS,Description="All filters passed">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (counting only informative reads out of the total reads) for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions for alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
    ##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
    ##FORMAT=<ID=GP,Number=G,Type=Float,Description="Phred-scaled posterior probabilities for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MB,Number=4,Type=Integer,Description="Per-sample component statistics to detect mate bias">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=PRI,Number=G,Type=Float,Description="Phred-scaled prior probabilities for genotypes">
    ##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias">
    ##FORMAT=<ID=SQ,Number=A,Type=Float,Description="Somatic quality">
    ##DRAGENCommandLine=<ID=HashTableBuild,Version="SW: 01.003.044.3.9.5, HashTableVersion: 8",CommandLineOptions="dragen --build-hash-table true --output-directory grch37 --ht-reference grch37.fasta">
    ##DRAGENCommandLine=<ID=dragen,Version="SW: 05.121.645.4.0.3, HW: 05.121.645",Date="Thu Mar 02 09:58:00 UTC 2023",CommandLineOptions="-f -c output.cfg -r /references/grch37">
    ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
    ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
    ##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
    ##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
    ##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (informative and non-informative); some reads may have been filtered based on mapq etc.">
    ##INFO=<ID=FractionInformativeReads,Number=1,Type=Float,Description="The fraction of informative reads out of the total reads">
    ##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
    ##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
    ##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
    ##FILTER=<ID=DRAGENSnpHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=DRAGENIndelHardQUAL,Description="Set if true:QUAL < 3">
    ##FILTER=<ID=LowDepth,Description="Set if true:DP <= 1">
    ##FILTER=<ID=PloidyConflict,Description="Genotype call from variant caller not consistent with chromosome ploidy">
    ##FILTER=<ID=base_quality,Description="Site filtered because median base quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=filtered_reads,Description="Site filtered because too large a fraction of reads have been filtered out">
    ##FILTER=<ID=fragment_length,Description="Site filtered because absolute difference between the median fragment length of alt reads and median fragment length of ref reads at this locus exceeds threshold">
    ##FILTER=<ID=low_af,Description="Allele frequency does not meet threshold">
    ##FILTER=<ID=low_depth,Description="Site filtered because the read depth is too low">
    ##FILTER=<ID=low_frac_info_reads,Description="Site filtered because the fraction of informative reads is below threshold">
    ##FILTER=<ID=low_normal_depth,Description="Site filtered because the normal sample read depth is too low">
    ##FILTER=<ID=long_indel,Description="Site filtered because the indel length is too long">
    ##FILTER=<ID=mapping_quality,Description="Site filtered because median mapping quality of alt reads at this locus does not meet threshold">
    ##FILTER=<ID=multiallelic,Description="Site filtered because more than two alt alleles pass tumor LOD">
    ##FILTER=<ID=non_homref_normal,Description="Site filtered because the normal sample genotype is not homozygous reference">
    ##FILTER=<ID=no_reliable_supporting_read,Description="Site filtered because no reliable supporting somatic read exists">
    ##FILTER=<ID=panel_of_normals,Description="Seen in at least one sample in the panel of normals vcf">
    ##FILTER=<ID=read_position,Description="Site filtered because median of distances between start/end of read and this locus is below threshold">
    ##FILTER=<ID=RMxNRepeatRegion,Description="Site filtered because all or part of the variant allele is a repeat of the reference">
    ##FILTER=<ID=str_contraction,Description="Site filtered due to suspected PCR error where the alt allele is one repeat unit less than the reference">
    ##FILTER=<ID=too_few_supporting_reads,Description="Site filtered because there are too few supporting reads in the tumor sample">
    ##FILTER=<ID=weak_evidence,Description="Somatic variant score does not meet threshold">
    ##ML=PackageVersion=4,ModelVersion=4,ModelCRC=2921457924
    ##contig=<ID=1,length=249250621>
    ##contig=<ID=2,length=243199373>
    ##contig=<ID=3,length=198022430>
    ##contig=<ID=4,length=191154276>
    ##contig=<ID=5,length=180915260>
    ##contig=<ID=6,length=171115067>
    ##contig=<ID=7,length=159138663>
    ##contig=<ID=8,length=146364022>
    ##contig=<ID=9,length=141213431>
    ##contig=<ID=10,length=135534747>
    ##contig=<ID=11,length=135006516>
    ##contig=<ID=12,length=133851895>
    ##contig=<ID=13,length=115169878>
    ##contig=<ID=14,length=107349540>
    ##contig=<ID=15,length=102531392>
    ##contig=<ID=16,length=90354753>
    ##contig=<ID=17,length=81195210>
    ##contig=<ID=18,length=78077248>
    ##contig=<ID=19,length=59128983>
    ##contig=<ID=20,length=63025520>
    ##contig=<ID=21,length=48129895>
    ##contig=<ID=22,length=51304566>
    ##contig=<ID=X,length=155270560>
    ##contig=<ID=Y,length=59373566>
    ##contig=<ID=MT,length=16569>
    ##contig=<ID=GL000207.1,length=4262>
    ##contig=<ID=GL000226.1,length=15008>
    ##contig=<ID=GL000229.1,length=19913>
    ##contig=<ID=GL000231.1,length=27386>
    ##contig=<ID=GL000210.1,length=27682>
    ##contig=<ID=GL000239.1,length=33824>
    ##contig=<ID=GL000235.1,length=34474>
    ##contig=<ID=GL000201.1,length=36148>
    ##contig=<ID=GL000247.1,length=36422>
    ##contig=<ID=GL000245.1,length=36651>
    ##contig=<ID=GL000197.1,length=37175>
    ##contig=<ID=GL000203.1,length=37498>
    ##contig=<ID=GL000246.1,length=38154>
    ##contig=<ID=GL000249.1,length=38502>
    ##contig=<ID=GL000196.1,length=38914>
    ##contig=<ID=GL000248.1,length=39786>
    ##contig=<ID=GL000244.1,length=39929>
    ##contig=<ID=GL000238.1,length=39939>
    ##contig=<ID=GL000202.1,length=40103>
    ##contig=<ID=GL000234.1,length=40531>
    ##contig=<ID=GL000232.1,length=40652>
    ##contig=<ID=GL000206.1,length=41001>
    ##contig=<ID=GL000240.1,length=41933>
    ##contig=<ID=GL000236.1,length=41934>
    ##contig=<ID=GL000241.1,length=42152>
    ##contig=<ID=GL000243.1,length=43341>
    ##contig=<ID=GL000242.1,length=43523>
    ##contig=<ID=GL000230.1,length=43691>
    ##contig=<ID=GL000237.1,length=45867>
    ##contig=<ID=GL000233.1,length=45941>
    ##contig=<ID=GL000204.1,length=81310>
    ##contig=<ID=GL000198.1,length=90085>
    ##contig=<ID=GL000208.1,length=92689>
    ##contig=<ID=GL000191.1,length=106433>
    ##contig=<ID=GL000227.1,length=128374>
    ##contig=<ID=GL000228.1,length=129120>
    ##contig=<ID=GL000214.1,length=137718>
    ##contig=<ID=GL000221.1,length=155397>
    ##contig=<ID=GL000209.1,length=159169>
    ##contig=<ID=GL000218.1,length=161147>
    ##contig=<ID=GL000220.1,length=161802>
    ##contig=<ID=GL000213.1,length=164239>
    ##contig=<ID=GL000211.1,length=166566>
    ##contig=<ID=GL000199.1,length=169874>
    ##contig=<ID=GL000217.1,length=172149>
    ##contig=<ID=GL000216.1,length=172294>
    ##contig=<ID=GL000215.1,length=172545>
    ##contig=<ID=GL000205.1,length=174588>
    ##contig=<ID=GL000219.1,length=179198>
    ##contig=<ID=GL000224.1,length=179693>
    ##contig=<ID=GL000223.1,length=180455>
    ##contig=<ID=GL000195.1,length=182896>
    ##contig=<ID=GL000212.1,length=186858>
    ##contig=<ID=GL000222.1,length=186861>
    ##contig=<ID=GL000200.1,length=187035>
    ##contig=<ID=GL000193.1,length=189789>
    ##contig=<ID=GL000194.1,length=191469>
    ##contig=<ID=GL000225.1,length=211173>
    ##contig=<ID=GL000192.1,length=547496>
    ##contig=<ID=NC_007605,length=171823>
    ##contig=<ID=hs37d5,length=35477943>
    ##reference=file:///references/grch37/reference.bin
    ##bcftools_filterVersion=1.9+htslib-1.9
    ##bcftools_filterCommand=filter -i 'TYPE="snp" ' /processing/sample.raw.vcf.gz; Date=Thu Mar  2 10:28:45 2023
    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  GFX0xxxxxx
    1   12783   .   G   A   6.16    PASS    AC=1;AF=0.5;AN=2;DP=37;FS=0;MQ=11.6;MQRankSum=-2.137;QD=0.45;ReadPosRankSum=2.878;SOR=1.949;FractionInformativeReads=0.973  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:10,26:0.722:36:4,14:6,12:4:38,0,1:6.1614,2.5231,7.0198:0,34.77,37.77:10,0,26,0:4,6,14,12
    1   13868   .   A   G   4.69    PASS    AC=1;AF=0.5;AN=2;DP=61;FS=2.275;MQ=16.39;MQRankSum=2.368;QD=0.36;ReadPosRankSum=-1.078;SOR=1.136;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:36,25:0.41:61:16,15:20,10:4:37,0,8:4.6918,2.1112,13.42:0,34.77,37.78:15,21,8,17:23,13,12,13
    1   13980   .   T   C   23.14   PASS    AC=2;AF=1;AN=2;DP=2;FS=0;MQ=5.7;MQRankSum=0;QD=6.5;ReadPosRankSum=0;SOR=2.303;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,2:1:2:0,1:0,1:19:61,24,0:23.143,20.887,0.056833:0,34.77,37.77:0,0,2,0:0,0,1,1
    1   14464   .   A   T   14.11   PASS    AC=1;AF=0.5;AN=2;DP=48;FS=7.243;MQ=22.38;MQRankSum=3.14;QD=0.86;ReadPosRankSum=2.527;SOR=0.267;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,37:0.771:48:8,16:3,21:7:77,6,0:14.113,0.98271,7.8591:0,34.77,37.77:3,8,19,18:4,7,20,17
    1   14907   .   A   G   3.67    PASS    AC=1;AF=0.5;AN=2;DP=77;FS=0.893;MQ=26.78;MQRankSum=-0.058;QD=0.41;ReadPosRankSum=5.956;SOR=0.532;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:48,29:0.377:77:31,13:17,16:4:36,0,22:3.6675,2.453,27.555:0,34.77,37.77:27,21,15,14:28,20,15,14
    1   14930   .   A   G   3.26    PASS    AC=1;AF=0.5;AN=2;DP=82;FS=3.051;MQ=27.04;MQRankSum=-0.363;QD=0.37;ReadPosRankSum=-0.277;SOR=0.754;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:51,31:0.378:82:30,16:21,15:3:35,0,21:3.2627,2.7897,26.606:0,34.77,37.77:28,23,14,17:29,22,17,14
    1   15211   .   T   G   6.66    PASS    AC=1;AF=0.5;AN=2;DP=70;FS=0;MQ=21.21;MQRankSum=0.533;QD=0.49;ReadPosRankSum=-0.905;SOR=0.743;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:28,42:0.6:70:12,22:16,20:6:40,0,11:6.6596,1.2235,15.27:0,34.77,37.77:15,13,23,19:15,13,17,25

Then a couple of chunks from random bits of the file:
offset 100000 lines:
    1   82279586    .   T   C   62.12   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.89;ReadPosRankSum=0;SOR=0.941;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,25:1:25:0,10:0,15:52:100,55,0:62.124,52.173,2.8992e-05:0,34.77,37.77:0,0,11,14:0,0,10,15
    1   82280153    .   G   A   62.73   PASS    AC=2;AF=1;AN=2;DP=31;FS=0;MQ=250;MQRankSum=0;QD=4.53;ReadPosRankSum=0;SOR=1.044;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,31:1:31:0,15:0,16:57:101,61,0:62.731,57.663,9.319e-06:0,34.77,37.77:0,0,18,13:0,0,14,17
    1   82280683    .   C   T   54.83   PASS    AC=2;AF=1;AN=2;DP=17;FS=0;MQ=250;MQRankSum=0;QD=5.77;ReadPosRankSum=0;SOR=1.061;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,17:1:17:0,8:0,9:45:93,49,0:54.832,45.914,0.00012581:0,34.77,37.77:0,0,10,7:0,0,11,6
    1   82280870    .   T   C   41.24   PASS    AC=1;AF=0.5;AN=2;DP=19;FS=1.808;MQ=250;MQRankSum=3.592;QD=2.58;ReadPosRankSum=1.858;SOR=1.14;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:11,8:0.421:19:8,2:3,6:40:76,0,45:41.236,0.00039866,47.817:0,34.77,37.77:4,7,2,6:5,6,6,2
    1   82280975    .   A   T   40.52   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.569;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=1.025;SOR=1.085;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:5,8:9,2:40:75,0,50:40.524,0.00040643,53.142:0,34.77,37.77:7,7,4,6:8,6,4,6
    1   82281492    .   T   C   44.14   PASS    AC=1;AF=0.5;AN=2;DP=22;FS=6.923;MQ=250;MQRankSum=3.873;QD=2.2;ReadPosRankSum=-0.334;SOR=2.064;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:13,9:0.409:22:9,4:4,5:44:79,0,51:44.143,0.00018328,54.3:0,34.77,37.77:6,7,7,2:6,7,6,3
    1   82282062    .   T   C   8.34    PASS    AC=1;AF=0.5;AN=2;DP=26;FS=1.974;MQ=247.84;MQRankSum=3.383;QD=1.04;ReadPosRankSum=-1.561;SOR=0.863;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:21,5:0.192:26:9,3:12,2:8:42,0,29:8.3401,0.69098,32.665:0,34.77,37.77:9,12,3,2:9,12,1,4
    1   82282378    .   C   T   42.61   PASS    AC=1;AF=0.5;AN=2;DP=24;FS=1.685;MQ=250;MQRankSum=4.07;QD=2.02;ReadPosRankSum=2.898;SOR=1.23;FractionInformativeReads=1  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,10:0.417:24:8,6:6,4:43:77,0,57:42.606,0.0002423,60.191:0,34.77,37.77:8,6,7,3:10,4,7,3
    1   82282788    .   C   T   48.99   PASS    AC=1;AF=0.5;AN=2;DP=30;FS=0;MQ=245.71;MQRankSum=2.889;QD=1.67;ReadPosRankSum=1.975;SOR=0.914;FractionInformativeReads=1 GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:14,16:0.533:30:6,9:8,7:46:84,0,46:48.985,0.00010665,49.217:0,34.77,37.77:8,6,10,6:6,8,7,9
    1   82282947    .   G   T   47.31   PASS    AC=1;AF=0.5;AN=2;DP=32;FS=10.606;MQ=250;MQRankSum=4.796;QD=1.56;ReadPosRankSum=0.68;SOR=0.725;FractionInformativeReads=1    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:15,17:0.531:32:6,8:9,9:46:82,0,48:47.312,0.00011286,51.287:0,34.77,37.77:11,4,7,10:8,7,8,9

offset 2000000 lines:
    8   142260995   .   T   C   65.76   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=248.93;MQRankSum=0;QD=4.69;ReadPosRankSum=0;SOR=1.179;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,28:1:28:0,14:0,14:53:104,56,0:65.759,52.946,2.3297e-05:0,34.77,37.77:0,0,11,17:0,0,13,15
    8   142261389   .   T   C   66.03   PASS    AC=2;AF=1;AN=2;DP=28;FS=0;MQ=250;MQRankSum=0;QD=4.67;ReadPosRankSum=0;SOR=0.997;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,28:1:28:0,12:0,16:54:104,57,0:66.026,54.004,1.812e-05:0,34.77,37.77:0,0,16,12:0,0,14,14
    8   142261541   .   G   A   65.34   PASS    AC=2;AF=1;AN=2;DP=30;FS=0;MQ=250;MQRankSum=0;QD=4.51;ReadPosRankSum=0;SOR=1.765;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,30:1:30:0,14:0,16:52:103,55,0:65.339,52.429,2.5886e-05:0,34.77,37.77:0,0,9,21:0,0,17,13
    8   142261618   .   A   G   64.46   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=250;MQRankSum=0;QD=4.57;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS   1|1:0,27:1:27:0,14:0,13:53:102,56,0:64.465,53.434,2.1227e-05:0,34.77,37.78:0,0,10,17:0,0,15,12:142261618
    8   142261621   .   C   T   60.33   PASS    AC=2;AF=1;AN=2;DP=25;FS=0;MQ=250;MQRankSum=0;QD=4.7;ReadPosRankSum=0;SOR=1.051;FractionInformativeReads=0.96    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB:PS 1|1:0,24:1:24:0,13:0,11:48:98,51,0:60.326,48.044,7.2481e-05:0,34.77,37.78:0,0,10,14:0,0,13,11:142261618
    8   142261967   .   C   T   68.75   PASS    AC=2;AF=1;AN=2;DP=37;FS=0;MQ=250;MQRankSum=0;QD=4.26;ReadPosRankSum=0;SOR=0.859;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,37:1:37:0,17:0,20:56:107,60,0:68.745,56.755,9.8367e-06:0,34.77,37.77:0,0,20,17:0,0,18,19
    8   142262327   .   A   G   68.48   PASS    AC=2;AF=1;AN=2;DP=42;FS=0;MQ=250;MQRankSum=0;QD=4.13;ReadPosRankSum=0;SOR=0.788;FractionInformativeReads=1GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB  1/1:0,42:1:42:0,21:0,21:56:106,60,0:68.476,56.786,9.8367e-06:0,34.77,37.77:0,0,22,20:0,0,11,31
    8   142262413   .   G   A   62.88   PASS    AC=2;AF=1;AN=2;DP=27;FS=0;MQ=231.86;MQRankSum=0;QD=4.75;ReadPosRankSum=0;SOR=1.302;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,27:1:27:0,16:0,11:51:101,54,0:62.881,51.183,3.5205e-05:0,34.77,37.77:0,0,10,17:0,0,9,18
    8   142263074   .   G   A   51.24   PASS    AC=2;AF=1;AN=2;DP=20;FS=0;MQ=143.27;MQRankSum=0;QD=5.31;ReadPosRankSum=0;SOR=0.892;FractionInformativeReads=1   GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,9:0,11:36:89,39,0:51.237,35.702,0.0012015:0,34.77,37.77:0,0,11,9:0,0,10,10
    8   142263171   .   C   T   52.9    PASS    AC=2;AF=1;AN=2;DP=22;FS=0;MQ=122.2;MQRankSum=0;QD=4.81;ReadPosRankSum=0;SOR=0.693;FractionInformativeReads=0.909    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    1/1:0,20:1:20:0,11:0,9:35:91,38,0:52.899,34.894,0.0014297:0,34.77,37.77:0,0,10,10:0,0,9,11

Download the dbSNP VCF file (latest version) and map to your input VCF using bcftools annotate. You'll want to match reference sequence version (hg19/hg38) and use the latest version of dbSNP to get the most up to date information.
dbSNP VCF FTP folder: https://ftp.ncbi.nih.gov/snp/latest_release/VCF/

Thanks, I eventually got this working although the brittleness with which it treats input file compression standards is rather annoying and took me ages to work out, hence my delay. I've finally annotated the vcf with rsid with bcftools, so thanks very much. Since my data is refed to GRCh37, and SNPedia is made with 38, although it works I wonder if it's worth 'getting' 38 somehow. Looks like generating it from the BAM or fastq would be the cleanest way to do this. Could I ask you to comment on that, and recommend any software or guides?

Yes, re-aligning to hg38 is the way to go. liftOver is an option too, but it has its shortcomings.

So regarding the realigning, is that a massive undertaking or relatively simple? I imagine it's simple - just testing the original files (bam or fastq etc) against the new reference genome hg38. Or does it require detailed info about the original sequencing operation that I am not likely to be able to get out of Dante (as in gringers reply re illumina snpchip versions being matched with specific dbsnp builds etc)?

Since it's WGS, my guess is that realigning should work just fine. Targeted sequencing might have created problems. You can ask @gringer for their take.

RamRS - thanks. Gringer, wonder if you come by these parts again if you could comment? As an extra ask, hoping its not too close to the wind to ask on SE about software recommendations, any suggestions for a tool to do this (I'm on Linux so Windows only software is out for me)?

As a partial answer to the specific questions about definitions:

rsIDs define a specific location and type of variation*, but the location of that variation can change within different reference genomes. This may include a change in chromosome, as well as a change in position.

Different dbSNP builds can add additional variation, but the same rsID will always refer to the same variant in all dbSNP builds that refer to that variation.

rsIDs can be invalidated in new dbSNP builds for various reasons. It's a good idea to always use the latest dbSNP build when annotating variants.

rsIDs are never recycled. There is a prescribed process for merging and splitting variants that ensures backwards compatibility. If a variation is merged, I think the approach is to always define the earliest-numbered variant as the primary variation. If a variation is split, I think the approach is to use a new rsID for each of the separate variations.

* by "variation", I'm referring to a type of modification in the DNA sequence, rather than a specific instance of that modification. This is explained in the dbSNP documentation (an "ss" is a submitted variant, a submission made to dbSNP):

Each ss that maps to the same position and is the same variant type is assigned to an existing RefSNP (rs) or assigned a new one.

Although that isn't entirely clear, this suggests that if a new variant is submitted that falls at the same position of an exiting one and is of the same type as the existing one, then the new one may simply be added to the same rsID as the existing one. So, for example, if we have a single nucleotide G => C variant at position N, and a new G => T variant is reported at the same position, then the new one could have the same rsID as the old one.

In response to an additional comment / question:

a post on the ANNOVAR site seems to suggest that no one really knows what rsid's are and that to different people they are different things, further muddying the waters.... Is this guy an authority on this or not?

I can't comment on authority for that guy; the authorship is unclear to me. What matters is how the dbSNP people understand dbSNP (i.e. the documentation link above), because they are the authority on it.
My understanding of it was helped because many years ago I attended a close-knit workshop that included a talk by Mike Feolo, and I had been working on a human dataset with many oddities, so I was able to ask lots of questions about corner cases. Here's the text from a section on rsIDs that I added into the appendix of my thesis; things have probably changed since then, but I expect the broad strokes are still correct:

Mike Feolo gave a talk at the [Working With The HapMap] Cambridge course on the dbSNP database. He talked about how new Single Nucleotide Polymorphisms (SNPs) were submitted, compared against currently existing SNPs, and merged if they were already in the database. The database RefSeq(rs) numbers are updated about twice a year, and whenever there is a new build of the human genome sequence. As a result of these builds, some rs numbers are merged to account for increased knowledge of mutations, a process that I was not aware of. This has an impact on my research, because I have received data from an Illumina 317k SNPchip, which I compare to the HapMap dataset. The rs numbers from the Illumina chip are for a specific dbSNP build, which may not be fully consistent with the most recent build (from which the HapMap data is retrieved). Mike has recommended that both the build number and the rs number are used when referring to a specific mutation. Sharma Buch pointed out that it was possible to ask Illumina for information about the merged rs numbers for their chips, which would allow datasets from different builds to be used together.

I think rsIDs are chr-pos specific, which is why multiple ALT alleles can exist in the same entry.

Thanks for all of this from all 3 of you. It's gradually helping me to form a picture. If anyone feels like reading and commenting on https://annovar.openbioinformatics.org/en/latest/articles/dbSNP/ it would be really useful. Is this guy an authority on this or not? His point that rsid means different thing to different people is revealing and if correct could prevent much head scratching and wasted hours reading and pondering this. It implies that one just has to get used to how rsid is used without trying to understand it too much (a sort of quantum mechanics like detachment...)

@Gringer your update is very interesting thanks.

Generally, consumer WGS companies do not deliver annotated VCF's. Annotations include rsIDs, gene names or similar area designations, and others.  But this can be done using the NIH dbSNP VCF files and the bcftools annotate command.  The annotation files at UCSC are modified for their GenomeBrowser and are often filtered and cleaned up.  While better and more consistent, they are usually not the latest release available from NIH.
Just make sure to match the build model and sequence (chromosome) naming convention between your annotation file and the original VCF. Dante usually delivers using the hs37d5 reference model. So use the .25 dbSNP VCFs files.  But dbSNP VCFs use accession names whereas Dante's hs37d5 reference model uses GRCh numeric ones.  So you will have to rename sequence names in one file or the other to make it work. Also, hopefully, the vendor has left aligned the InDels as that is what the dbSNP file has them defined as and what bcftools expects as well. Dante traditionally uses the DRAGEN pipeline with their Illumina sequencers; which does this. But they have recently switched to MGI DNB sequencers and the megaBOLT pipeline for customers outside the USA. We are not yet sure how they are processing InDels. Finally, Dante has historically delivered a separate SNP and InDel file.  You will often want to merge them or annotate them both separately then merge. bcftools can do this.
The dbSNP VCF file has well over 1 billion rsIDs now defined. And many fields of additional information that can be annotated in.  Often, if you just want to add a simple subset of the data, you can process the dbSNP VCF file to a simpler tab-separated (TSV, tab) file of entries. It can cut down the annotation file from 25 GB to 5 GB or less. See the bcftools annotation page for hints on using that tool to do this as well.
Personally, I find it easier to use online tools or services that have been tweaked.  Especially if they have nice browsers.  Sequencing dot com, usegalaxy.org and iobio.io are a few that come to mind.  Not all allow the export of the annotated VCF.
This and similar topics of utilizing your results are heavily discussed in the consumerWGS Facebook group; if you also wish to participate there.


text:
From SMILE get Amber force in Python
I have a list of SMILES of small molecules and I want to be able to simulate these molecules with an Amber force field in Python. Currently, I use RDkit to convert the smile into a PDB file:
from rdkit import Chem
from rdkit.Chem import AllChem

molecule = Chem.MolFromSMiles('smile')
molecule = Chem.AddHs(molecule)
AllChem.EmbedMolecule(molecule)
Chem.MolToPDBFile(molecule, 'molecule.pdb')

Then I use the ParmEd package:
import parmed as pmd

pdb = pmd.load_file('molecule.pdb')

This pdb object for 'CC(C)C' for example, looks like:
<Structure 14 atoms; 1 residues; 13 bonds; NOT parameterized>

Then I want to create a prmtop file, which is needed to initialise an Amber force:
pdb.save('molecule.prmtop', overwrite=True)

This throws the error:
RuntimeError: atom types are not assigned

I'm not sure why this error gets thrown, since the atom types do seem to be inside the PDB file. I'm new to trying to do these kinds of things, so I'm not completely sure if more information is needed.
If this method in general is not preferable, what is the best way to achieve this (preferably on python)?
Any help would be greatly appreciated.


text:
To find genes that don't change in RNA-seq, Deseq2 has altHypothesis="lessAbs". Is there a way to make limma do the same thing for proteomics?
I love that Deseq2 has altHypothesis="lessAbs" !!!!!!! I've used it a ton for RNA-seq. However, now I'm working with mass spect data using limma. Is there a way to make limma do the same thing?

Hi and welcome to the site! Could you please [edit] your question and clarify a bit? First, please make the title be a title: a short, concise description of what the question is about. Second, please give us a bit more context. I think limma is an R package, right? And what exactly would this option do? Identify genes that don't change? Between what and what? Are you comparing two sets of mass-spec data? Please add all extra info to your question, don't answer in the comments, just add it to the main question.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

There is no direct support for this by best knowledge. Similar questions have been asked previously in the Bioconductor support forum, for example:
https://support.bioconductor.org/p/64300/#64793
They discuss options around the TOST concept.
For edgeR:
https://support.bioconductor.org/p/66283/#66287
From these posts please go back and forth to find other posts in that forum to capture what the developers of limma/edgeR responded on the matter.
Depending on your exact purpose it might even be enough to go for genes with large AveExpr and a nominal (not adj) p-value close to one.


text:
Is my reference sequence too small?
I'm trying to map ONT long reads to a portion of a gene I'm looking at. The region is about 25bp long. When I search for the region in the document it pulls up the sequence in every read but when I map to reference with minimap2 it produces no contigs.
Is 25bp too small to be a reference? I can potentially extend it but it's right in an area of conservation in this region of a gene and I would rather not extend much farther.

Can you explain why you would want such a small sequence? If you're mapping, why not use the entire genome? I don't know the specific step size used by minimap2 but I would be very surprised if 25bp is enough, yes. Can you explain your use case in a bit more detail? Why would you want to use a reference that is shorter than even short reads, let alone long reads? What is the perceived benefit? Why is it that you would rather not extend?

If you don't have a reference genome for your organism/target, you should assemble the reads instead of trying to align them to a short conserved area. I think de novo assemblers exist that are developed specifically for Nanopore data.

That's not the right way to do it.
You either have reference or not. If not, you need to assemble one (that problem is called genome assembly). If yes, you map all the reads to the whole reference - the process of mapping is made for this - placing individual WGS reads to their corresponding locations and as a matter of fact, missing portions of the reference can mess it up (that was one of the reasons why there was a push for T2T for human).
Once the reads are mapped, you can use any genome browser (e.g. IGV) to zoom in to the exact location you are interested in and check the stack of reads mapping in there. Given how small the region is, manual inspection of the quality of the mapping is the best you got.


text:
file order after nextflow groupTuple
Below is the code from nextflow site (https://www.nextflow.io/docs/edge/operator.html#grouptuple).
For the output, the first is A B C, and the second is C, A.
Is there a way to make the output like this
[1, [A, B, C]]
[2, [A, C]]
That is, the stuff inside [] is in some order.
Channel
    .of(
        [1, 'A'],
        [1, 'B'],
        [2, 'C'],
        [3, 'B'],
        [1, 'C'],
        [2, 'A'],
        [3, 'D']
    )
    .groupTuple()
    .view()

[1, [A, B, C]]
[2, [C, A]]
[3, [B, D]]

You can change the sorting criteria for the grouped items using the sort keyword. To perform a natural sort, just use sort: true. For example:
Channel.of(
    [1, 'A'],
    [1, 'B'],
    [2, 'C'],
    [3, 'B'],
    [1, 'C'],
    [2, 'A'],
    [3, 'D']
)
.groupTuple(sort: true)
.view()

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.04.0
Launching `main.nf` [special_darwin] DSL2 - revision: 2e641ccfd9
[1, [A, B, C]]
[2, [A, C]]
[3, [B, D]]

Otherwise, a custom sorting criteria can be specified using a closure.


text:
Calculating total number of Unique gene in every genome from ROARY pangenome output
I performed a ROARY from 100 genome. I got several file including gene_presence_absence.csv  from the results. I want to know how many unique gene(unshared gene, based on blast homology, same gene (duplicate gene) may be unique if no homology is there between them) from presence absence csv file or other file. Kindly suggest me how to calculate total number of unique gene for each genome which is only present in one genome

How do you define "unique genes"? Unique to what? You mean genetic sequences that code for protein which are only present once in a given genome? Do you mean unique to that species? Unique how? In terms of sequence? In terms of function? How much similarity do you allow? What species are you working with? Please [edit] your question and explain what you need and also what kind of data you have more clearly.

Thank you for the reply, based on presence absence  gene data file, only gene found in one genome but absent in other genome is considered as unique gene, ROARY use use presence absence data based on blast homology of protein within selected protein so, at end I want to get only unique gene from the table generated by PROKKA in gene presence absence csv file.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Once I did ROARY, I got several files, right, there will be some shared gene, core gene , and unique gene, right, I want to know how many unshared gene (unique gene) are present in every genome.

Please [edit] your question and provide us with an example input file and the output you would expect.

I am trying to understand ROARY output results , sure I will do

Did you try query_pan_genome?
Ref. https://sanger-pathogens.github.io/Roary/

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Given a FASTA file with a set of genes, how do I determine that set's conservation among a genetic order?
Given a set of genes in a FASTA-formatted file, for example: exFasta.fa, what is the best way to determine how well conserved the genes in that file are among an order of organisms, for example: Orthoptera?
BLAST is an obvious option, but I am not sure that it is the most efficient way considering there are likely thousands of organisms in that order. OrthoFinder also comes across as a valid option, but I am not sure how to implement such a thing considering the required inputs.
Any advice would be greatly appreciated!

Are these annotated, known genes? Or do you also need to detect/identify the potential homologs in the target species?

Orthoptera (grasshoppers, crickets, locust) will be annotated (mostly). What they want to do is skim the draft genomes of Orthoptera. I like @gringers answer, it's a good one. As a project it would work.

It's possible that MOLE-BLAST could help for you:

MOLE-BLAST is an experimental tool that helps taxonomists find closest database neighbors of submitted query sequences. It computes a multiple sequence alignment (MSA) between the query sequences along with their top BLAST database hits, and generates a phylogenetic tree. Query sequences in the tree are denoted with highlighted node labels. If the input sequences come from different genes or loci, MOLE-BLAST can cluster them and compute an MSA and a phylogenetic tree for each locus separately. This tool may eventually be incorporated into our supported suite of BLAST applications based on the feedback we receive.

I haven't actually used MOLE-BLAST, but I have used Smart BLAST, which does a somewhat similar phylogenetic search (although on a restricted diverse taxonomic set) for protein queries.

thank you! I'll try that.


text:
PROKKA setupdb / cmpress error
I am trying to run --setupdb in Prokka, but I'm getting this error that I can't seem to find the answer for.
When I run
prokka --setupdb
I get
[18:57:31] Determined cmpress version is 001001 from '# INFERNAL 1.1.4 (Dec 2020)'
[18:57:31] Pressing CM database: /home/carol/.local/lib/prokka/db/cm/Viruses
[18:57:31] Running: cmpress \/home\/carol\/\.local\/lib\/prokka\/db\/cm\/Viruses

Error: File existence/permissions problem in trying to open CM file /home/carol/.local/lib/prokka/db/cm/Viruses.
CM file /home/carol/.local/lib/prokka/db/cm/Viruses not found (nor an .i1m binary of it)

[18:57:31] Could not run command: cmpress \/home\/carol\/\.local\/lib\/prokka\/db\/cm\/Viruses

Does anyone know how I can fix it?

This is a known problem with Prokka and was previously answered here https://bioinformatics.stackexchange.com/questions/20620/barrnap-bacterial-rrna-predictor-script-permission-denied-error-when-running-pro/20621#20621
Please let us know if this is okay, then we can close the question.


text:
Cut and Tag analysis
I am trying to do Cut and Tag analysis for a set of files and I want to perform a procedure called Tn5 normalization of samples.
The basic goal is to generate bedgraph files and does normalization with the weights of a sequencing depth. I have a set of function files that I have defined and then a set of Tn5 normalization scripts. The scripts are both quite long and hence I will be sharing parts of the scripts to get hold of my error and also provide as much information as possible.
My bedgraph files look like this:
replicate species Msample serotype protein            sample      condition    group
1        R1     NHP      02     LK03      NA R1-NHP-02-LK03-NA R1-NHP-LK03-NA LK03-NHP
2        R1     NHP      04     LK03      NA R1-NHP-04-LK03-NA R1-NHP-LK03-NA LK03-NHP
3        R1     NHP      13      KP1      NA R1-NHP-13-KP1-NA R1-NHP-KP1-NA KP1-NHP
4        R1     NHP      42      KP1      NA R1-NHP-42-KP1-NA R1-NHP-KP1-NA KP1-NHP
5        R1     NHP      55     LK03      NA R1-NHP-55-LK03-NA R1-NHP-LK03-NA LK03-NHP
6        R1     NHP      57      KP1      NA R1-NHP-57-KP1-NA R1-NHP-KP1-NA KP1-NHP

I must load a complete set of files called RDS which basically bears a list of the individual bed graph files.
The read lengths for each of the files look like this:
seqnames start  end normcount
73694699     opt2     0   52 0.0000000
73694700     opt2    52   56 0.0243838
73694701     opt2    56   83 0.0731515

The Tn5 R normalization scripts which are relevant
library(dplyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(gridExtra, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(Gviz, lib.loc = "/home/ag111/R/x86_64-pc-linux-gnu-library/4.1") #needed to be installed
library(ggplot2)
library(RcppArmadillo, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
#library(stringi, lib.loc = "/scg/apps/software/r/4.0.3/lib64/R/library")
library(tidyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(RColorBrewer, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(ggpubr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library") #needed to be installed
library(GenomicFeatures, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(GenomicRanges, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")

source('AG_CnRT_functions.r')
dir.create(file.path("Tn5exITR_InVivo_norm_full"))

intermediates_dirs <- c(SeqDepth='intermediates/')
#intermediates_dirs <- c(SeqDepth='/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/')

normalization_method <- 'SeqDepth'

AAV_gsize <- 4410
AAV_anno <- read.delim('AAV_opt2_anno.txt')

#figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

############## Tn5 ###########
AAV_normbedgraphs <- readRDS("/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/AAV_rmdup_SeqDepthNormbedgraphs.rds")

naming_order <- c('replicate','species',"Msample",'serotype','protein')

AAV_normbedgraphs_exITR <- AAV_normbedgraphs %% lapply(function(x){
  df <- x %% subset(start = 185 & end <=  4167)
  return(df)
})

AAV_normcov_regions_meta_exITR <- list()
AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG(
  AAV_normbedgraphs=AAV_normbedgraphs_exITR,
  naming_order=naming_order,
  AAV_gsize=AAV_gsize,
  AAV_anno=AAV_anno)

Merging data
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "AM"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("AM", "LK03", "DJ"))

AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "KP1"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("KP1", "LK03"))

metadata <- get_metadata_AG3(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             #naming_order=naming_order)
metadata <- get_metadata_KP(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             naming_order=naming_order)

head(metadata)

#coverage separated by AAV regions and species
#to_plot <- c()
#to_plot <- AAV_normcov_regions_meta_exITR[[normalization_method]] %% gather('region','normcov',2:6)
#to_plot$region <- factor(to_plot$region, levels=c(as.character(AAV_anno_exITR$region),'AAV_normcov'))
#max_y <- ceiling(max( log2(to_plot$normcov1) ) )
#min_y <- floor( min( log2(to_plot$normcov1) ) ).

After loading of the rds file the error I am getting while I am executing the code block
 AAV_normcov_regions_meta_exITR <- list()
 #AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG3(
 AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG(
   AAV_normbedgraphs=AAV_normbedgraphs_exITR,
   naming_order=naming_order,
   AAV_gsize=AAV_gsize,
   AAV_anno=AAV_anno)

Error in if (!all(AAV_normcov_regions$sample == metadata$sample)) { :
  missing value where TRUE/FALSE needed.

The code block which I think for the function file that is valid  is
AAV_normcov_regions[is.na(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.null(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.na(AAV_normcov_regions)] = 0
  AAV_normcov_regions <- AAV_normcov_regions[ match(metadata$sample,AAV_normcov_regions$sample) , ]
  if ( !all(AAV_normcov_regions$sample == metadata$sample) ){
    stop('Unsuccessful matching to metadata')
  }
  return(AAV_normcov_regions).

This function is working for another set of samples to generate a plot list.
The plot list I am having is to_plot
[1] sample      AAV_normcov replicate   species     Msample     serotype    protein     condition   group
<0 rows>

(or 0-length row.names) where as what is expected is something like this
  sample   Enhancer        FLuc      polyA   Promoter       WPRE AAV_normcov replicate batch species drug serotype condition group
1  R1-LiD-H-N-DJ 0.42244801 0.290929537 0.43640203 0.22379898 0.34125606  0.26650388        R1   LiD       H    N       DJ      H-DJ  H-DJ
2  R1-LiD-H-N-iG 1.42930414 0.259549496 0.57248120 0.37425687 0.42189340  0.36092640        R1   LiD       H    N       AM      H-iG  H-iG
3  R1-LiD-H-N-LK 0.25754367 0.317786560 0.27373645 0.19470267 0.29818285  0.24753280        R1   LiD       H    N     LK03      H-LK  H-LK
4  R1-LiD-M-N-DJ 0.24578258 0.097169668 0.14489638 0.06470209 0.11091813  0.09308374        R1   LiD       M    N       DJ      M-DJ  M-DJ

I am happy to share other code blocks if it is required to answer the question.

Please see our help page on formatting. I tried to format your code as code, but please check and make sure I didn't miss something.

Is there a typo in your code? You have a x %% subset which should probably be x %>% subset. Did the code run or did you edit it after copy-pasting and make a mistake there?

It will be indeed be x %>% subset. There was a mistake during copy pasting.

Check your code for more copy-paste errors and make sure everything lines up. Edit your post and fix all the errors you find, including the %% subset error. Please use proper code formatting for inline code and code blocks.

The file that you include as a bedgraph example doesn't look like bedgraphs.


text:
STACKS: process_radtags only reads 1 input file and returns >90% barcode not found drops
I am trying to demultiplex some paired-end ddradseq data and am running into an issue with STACKS in that the program only seems to read 1 of my files for input (there are 2) and also results in over 90% of my barcodes identified as not found, and therefore dropped. As a result, my retained reads are extremely low (1.1%). I ran this same script before on the same dataset and had good results so I am not sure why I am coming across this issue now?
The current script I am using is as follows:
module load stacks

process_radtags -p /home/heather/scratch/raw_inputs/ -P -o /home/heather/scratch/processrad_out \
 -b /home/heather/scratch/barcodes/lesp_barcodes_batch1.txt \
 -c -q -r --inline_inline --renz_1 sbfI --renz_2 mluCI

echo "Process radtags complete. Check log file to see dropped reads then proceed to alignment"

My input files are named HFL_R1_001.fastq.gz and HFL_R1_002.fastq.gz respectively (I have simplified these for the sake of this post, but the correct file names are specified in my scripts so this is not the issue). The program runs, but this is the output of the log file that I get which shows very low retained reads.
BEGIN per_file_raw_read_counts
File    Retained Reads  Low Quality     Barcode Not Found   RAD cutsite Not Found   Total
HFL_S1_L001_R1_001.fastq.gz 11092204        3150    1029697796  12760620        1053553770
END per_file_raw_read_counts

BEGIN total_raw_read_counts
Total Sequences 1053553770
Barcode Not Found   1029697796  97.7%
Low Quality     3150    0.0%
RAD Cutsite Not Found   12760620        1.2%
Retained Reads  11092204        1.1%
END total_raw_read_counts

I have also tried running a script where I specified with the -1 and -2 flags which files I wanted the program to input specifically and get the same result, with the difference of the program seems to read data from both files, but still only finds the 1 input file:
Processing paired-end data.
Using Phred+33 encoding for quality scores.
Found 1 paired input file(s).
Searching for single and paired-end, inlined barcodes.
Loaded 191 barcodes (4-6bp / 6bp).
Will attempt to recover barcodes with at most 1 / 1 mismatches.
Processing file 1 of 1 [HFL_S1_L001_R1_001.fastq.gz]
  Reading data from:
  /home/heather/scratch/raw_inputs/HFL_S1_L001_R1_001.fastq.gz and
  /home/heather/scratch/raw_inputs/HFL_S1_L001_R2_001.fastq.gz

The rest of the log file from this particular script looks almost identical to the one described above, where retained reads are around 1%.
I have checked my barcode file and there are no unreadable characters and I can't imagine why I would be getting so many lost barcodes when this ran previously with no problem? If anyone has any insight it would be much appreciated. I am relatively new to bioinformatics and don't really know where to go from here.