cannin/biostars_qa huggingface dataset:
Question: In the [VCF format][1], there is the option for the ID field to have multiple semi-colon separated values. In theory, there could be two dbSNP rsIDs in a single line (i.e. two indels at chr:pos), but for programming purposes, that should not happen, correct? dbSNP has merged all variants for a given position to a common rsID?

[1]: https://samtools.github.io/hts-specs/VCFv4.3.pdf

Answer: > dbSNP has merged all variants for a given position to a common rsID?

I'm afraid no:

$ wget -q -O - "ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/All_20180418.vcf.gz" | gunzip -c | grep -v "#" | cut -f 1,2 | uniq -d | head
1	10051
1	10055
1	10108
1	10109
1	10128
1	10132
1	10177
1	10228
1	10229
1	10235
.

$ wget -q -O - "ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/All_20180418.vcf.gz" | gunzip -c | grep -v "#" | cut -f 1,2,3,4,5 | grep -w 10051 -m2
1	10051	rs1052373574	A	G
1	10051	rs1326880612	A	AC
Question: Hello all,

I was wondering if there is a guide to modern genome assembly - that is, what often involves both short and long reads, as well as Hi-C sometimes?
What tricks can be used for typical problems (contaminated reads, highly heterozygous genome, etc)?

Thank you in advance, as always

Answer: There are some guidelines from the Earth BioGenome Project here: https://www.earthbiogenome.org/assembly-standards

These address quality requirements and standards rather than how to achieve them, but they also mention the [Vertebrate Genome Project assembly pipeline][1]. It might be useful to look at the tools VGP pipeline use and related publications, e.g. the [Rhie et al., 2021 paper][2]. Indeed, the Methods and Supplementary might sufficiently address the issues you mention.

What helps against highly heterozygous genome, btw? High levels of inbreeding over many, many generations (but that has maybe fallen out of fashion).

[1]: https://github.com/VGP/vgp-assembly
[2]: https://www.nature.com/articles/s41586-021-03451-0
Question: I have a tree with more than 1000 sequences. I want to colour branches based on a string in the tip labels.

For example, all the branches leading to tip labels that have the word 'bacteria' in, I want to label blue, and those with the word 'archaea' I want to label red.

It is straight forward to do this for the tip labels themselves, just by using ctrl + f, and searching for your string of interest, then selecting a colour, but it seems like it is not possible to also colour the branches. Does anybody know a straightforward way of doing this, ideally without needing to write a script?

Cheers!

Answer: In FigTree you can just click on the nodes or branches of interest and colour them accordingly using the colour wheel.

If you change your selection mode to Clade, all child nodes and branches will be selected at the same time, so you only need to select the common ancestor for any set of branches you want to colour.

I don't believe there is any way to do this 'automatically' based on a given string, but it is fairly trivial to do by hand.

Example:

1. Change selection mode to Clade (Blue box)

2. Select a node/branch set (black box)

3. Apply a colour (red box)

![enter image description here][1]

[1]: https://s8.postimg.cc/6k30b5stx/Screen_Shot_2018-07-31_at_16.15.45.png
Question: I'm a Computer Science undergraduate writing my thesis which is about analyzing open source bioinformatics projects and extracting standalone modules. I need those projects to be written in Java and preferably be stable or mature. So far I've found a bunch of projects via SourceForge and I was wondering if there are any notable ones that are not hosted there. A different repository that has a bioinformatics category would be helpful too.

Here's a list of the most promising projects I have so far:

- Jenetics
- JGap
- Jmol
- Juicebox
- Picard
- RDP Classifier
- Shim
- Tassel
- The Chemistry Development Kit
- Toxfree
- TreeView
- VarScan

Thanks in advance.

Answer: [JEnsembl][1] is a JAVA API for Ensembl.

[Contrail][2] and [Cloudbase][3] are assemblers. The latter got some new life when an enterprising programmer [retooled it for GPUs][4].

[1]: http://jensembl.sourceforge.net/
[2]: http://contrail-bio.svn.sourceforge.net/
[3]: https://github.com/ice91/CloudBrush
[4]: https://github.com/zenlc2000/CloudbrushGPU
Question: Hi, guys! I am a newbie in pabio sequence analysis , recently I download some pacbio RSII data from ncbi sra database, after fastq-dump , I get a fastq file which I think is the raw subread data.

But now I want to do ccs(Circular Consensus Sequence calling) step , the pacbio ccs binary software need the bam or bax (bax can be transformed to bam )file as input , also in the sra website , there is no h5 format file submitted by submitter , so there only has fastq file , so I don't know how to do with the pacbio fastq file ? Thanks!

Answer: Since you say that you are getting Fastq/fasta files, the first step is to correct them. The best way to go about it is to either use Falcon till the generation of preads or Canu's correct-trim mode or use proovread.

If you had the bax/bam files, this would have been easier as you could import them into SMRTlink and you could have used it's pushbutton nature :)
Question: Hello,

I am thinking of designing an algorithm to rank variants in order of most to least likely to have a functional effect in a autoimmune disease.

The problem is, if I wanted to train such an algorithm, it would be helpful to have a list of variants that are causative, or almost certainly causative, of an autoimmune disease.

For instance, in the Rheumatoid Arthritis literature, there are a lot of associated SNPs, but we are reasonably sure that amino acid positions 71 and 74 of HLA-DRB1 actually CAUSE the disease, rather than being an associated SNP in LD with a SNP that actually increases or decreases risk of the disease.

It can be any autoimmune disease, but mendelian diseases etc will not help. The goal will be to mine characteristics of these causative SNPs.

I am aware of algorithms like CADD, PICS, etc., I am more in need of a list of validated causal snps than anything else.

Does anyone know of such a source?

-------------ADDENDUM - someone posted below the NHGRI catalog-------------------

The problem is that the NHGRI catalog lists associations, not (necessarily) causal variants.

To determine a causal variant, we need follow up experiments in the wetlab to show necessity and sufficiency. It is thought that the lead SNP is actually causal in 10% or fewer of cases (Farh et al 2015).

I am looking for a resource only for the latter type of SNP.

Answer: Maybe https://opensnp.org or http://www.snpedia.com/ could help you.
Question: Hi, I have 10 microarray data sets (each data set related to a disease) in which I already compiled to get differentially expressed genes (DEGs) in every individual data set. I want to extract common DEGs of all the datasets, is there any tool/R package/R function to do it.

I have the ProbesetIDs, Genenames, GeneSymbols, Entrez IDs, LogFC values, B values, P values; in all the datasets.

Any help would be appreciated. Thanks!

Answer: Trying to do something like this?

```
gene <- c("a", "b", "c"); expr <- c(2, 2, 3);
x <- data.frame(gene, expr)
gene <- c("c", "b", "e");
y <- data.frame(gene, expr)
temp <- merge(x,y,by=match("gene", colnames(x)))
gene <- c("c", "e", "d");
z <-data.frame(gene, expr)
temp <- merge(temp, z, by=match("gene", colnames(temp)))
#...
#repeat for all your dataframes
```

Essentially you're matching the dataframes 1 by 1 on a specific column until you've merged all of them.

In your case you could go matching by e.g. "Genenames".
Question: Dear all,

I have the following experimental situation: a series of normal samples transcriptomes, and a series of aberrant transcriptomes.
However, I don't want to perform the standard analysis Test Vs Normal.
Instead, I would like to be able to call for each aberrant sample which genes are deregulated and which genes are not.
I know this is fragile from a statistical point of view, having no replicas for my aberrant samples.
Can I use standard tools for DEG analysis or is there a strategy particularly amenable to achieve this result?

The aberrant samples are cancers and for the scope of my investigation, I don't want to assume that they are equal.
Even if that means losing some power, I want to look at sample specific deregulation comparing that with the distribution in normal samples.

I hope my question is not trivial since I searched for a while and I didn't find a proper answer.
Thank you in advance,

Gian

Answer: You can use standard DE tools as long as you have replicates for one of the groups. The power is not going to be great, you'll need more controls, and it's questionable how biologically meaningful the results really are, but it'll work.
Question: Hi everybody,

I am running GATK pipeline for variant calling in trinity with following command.

/usr/local/bin/Trinityrnaseq-v2.6.6/Analysis/SuperTranscripts/AllelicVariants/run_variant_calling.py --st_fa ./SuperDuper.fasta --st_gtf ./SuperDuper.gff -p ./FCHG1.fq.gz ./FCHG.fq.gz -o ./variant_calls_outdir

But I got following error:

Error, missing path to Picard-Tools in $PICARD_HOME.
based on the instructions on the trinity web, I need to set the installation directory to the environmental variable ${PICARD_HOME}. for this purpose (adding the picard path to $PATH) I tried following command but my problem was not solved.

PATH="$PATH:$HOME/picard-2.23.0"
echo $PATH
output

/home/ubuntu/bin:/home/ubuntu/.local/bin:/home/ubuntu/anaconda3/bin:/home/ubuntu/perl5/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/home/ubuntu/picard-2.23.0

I don’t know how can I set the installation directory to the environmental variable ${PICARD_HOME}.

best regards

Answer: Please read the error message properly. You need to set `PICARD_HOME` to the installation directory, not add it to `PATH`. Try

PICARD_HOME=$HOME/picard-2.23.0

and then run the script again.
Question: I have a map file and gff file and all I want to do is for each line in the map file match the first column for each line in gff file and any matches and replaced with the second column from map file. I've tried the map id script from maker2 but isn't working beyond the mRNA field, the parent CDS field is not changed so not much use to me.

gff file tab delineated:

```
Chromosome_1 Geneious mRNA 421 2049 . - . ID=12;Parent=12-gene
Chromosome_1 Geneious gene 421 2049 . - . ID=12-gene
Chromosome_1 Geneious CDS 421 2049 . - . ID=12:cds;Parent=12
Chromosome_1 Geneious CDS 3747 6569 . - . ID=rmRNA-1:cds;Parent=rmRNA-1
Chromosome_1 maker mRNA 3747 6569 . - . ID=rmRNA-1;Parent=gene1
Chromosome_1 maker gene 3747 6569 . - . ID=gene1
```

Map file tab delineated

```
12 B1
rmRNA-1 B2
```

output desired:

```
Chromosome_1 Geneious mRNA 421 2049 . - . ID=B1;Parent=B1-gene
Chromosome_1 Geneious gene 421 2049 . - . ID=B1-gene
Chromosome_1 Geneious CDS 421 2049 . - . ID=B1:cds;Parent=B1
Chromosome_1 Geneious CDS 3747 6569 . - . ID=B2:cds;Parent=B2
Chromosome_1 maker mRNA 3747 6569 . - . ID=B2;Parent=gene1
Chromosome_1 maker gene 3747 6569 . - . ID=gene1
```

Answer: This is a little python script that will do the job.

```py
gff_file = 'gff.txt'
map_file = 'map.txt'
out_file = 'output.txt'

my_map = dict()
with open(map_file, 'r') as f:
for line in f:
data = line.split('\t')
my_map.update({data[0]:data[1]})

with open(out_file, 'w') as out:
with open(gff_file, 'r') as f:

for line in f:

data = line.split('\t')

for key, value in my_map.items():
if key in data[-1]:
print(key,value)
print(data[-1])
data[-1] = data[-1].replace(key, value.strip())
print(data[-1])

for item in data:
out.write(item+'\t')
```
Question: Hi All,

Apologies for the technical/irrelevant question.

Some gene names start with APR*/MARC*/SEPT* etc default converted into date format.

How to avoid this automatics conversion of gene symbols to date in Excel.

Thanks..

Answer: If we assume that Excel is an unavoidable evil (and for most this is unfortunately a credible assumption)...

When importing your delimited data into Excel, be sure to set the column type appropriately. The default for all columns is 'General' - this will lead to the inevitable auto-correct mangling. Gene name columns should be set to 'Text', then the septins etc will be left unmolested.

Once mangling has occurred, correction needs to be approached carefully - the main issue with what Excel does here is that it *changes the underlying data*, rather than only changing things in the display layer. This means your gene names are turned into floats (which is Excel's underlying representation of dates). Of course, Excel being Excel, the float used is different depending on the platform: https://support.microsoft.com/en-us/kb/214330 (also some versions of Excel think 1900 was a leap year, others know better - this is a very deep rabbit hole...). You could conceivably write some code to handle the conversion of dates back to gene names, but it would have to consider these caveats and edge cases veeeerrrry carefully.

All in all, I agree with Michael (although accept that this is not always possible).
Question: I saw some research showed that blast is 15% more sensitive than diamond, while I am interested in using blastp in OrthoFinder, I can't seemed to find the parameter of Blast being documented, there is also nothing in the log.file from Orthofinder output. How did blastp from Orthofinder select an ortholog? based on which parameter cutoff? I believe using OrthoFinder is better than bidirectional blast in finding orthologs, but I hope I can find more documentation that support this thought. The only input parameter Orthofinder accepted is -S for selecting propgram of interest( is there any other blast option ? ) . Any advice, suggestion, post that I have missed from biostar or criticism are welcome. Thank you in advance.

Answer: Just to put the [other answer's][1] claim in context: the latest version of `DIAMOND` (`v2.0.7`) is as sensitive as `blastp` while being (at least) 80x faster. I believe `OrthoFinder` also supports `MMseqs2`, and there you'd get a 10-12x speedup at `blast`-like levels of sensitivity. `OrthoFinder v2.5.2` uses `DIAMOND v2.0.6` but that should still be somewhere in this ballpark. So there's no point NOT going with `DIAMOND` or `MMseqs2` unless you want everything to be really slow for some reason.

`OrthoFinder` first takes all matches with an **e-value of 10^-3** (not sure if this is still the cut-off) from an all-vs-all search of the input sequence sets. Then it carries out a normalization procedure that ensures that all the best hits between all the input species have the same scores irrespective of sequence length and/or phylogenetic distance. Thereafter, for each query gene, the "worst" reciprocal best hit/`RBH` (or Reciprocal Best length-Normalised hit/`RBNH`) is identified. All matches that have a score better than this match are taken to form the orthogroup for that query, irrespective of whether they are `RBNH`s or not. These sets of query-matches are then clustered using the MCL algorithm to yield the orthogroups. Then gene trees are constructed and resolved using these orthogroups, leading to the identification of orthologs (and paralogs, albeit indirectly). More details can be found [here][2], [here][3], and [here][4]. (I have to admit that the `OrthoFinder` documentation was hard to follow, and I actually had to spend quite a bit of time in the `GitHub issues` to figure out a couple of things for myself.)

So in a sense, yes, `OrthoFinder` is better than just reciprocal best hits/`RBH`s (e.g., bi-directional `BLAST`) because the latter does not perform all these additional procedures. And it uses the 10^-3 e-value cut-off to start off with a broad set of matches that can then have the false positives sieved out of them in the subsequent steps.

[1]: https://www.biostars.org/p/9475302/#9475310
[2]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0721-2
[3]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1832-y
[4]: https://github.com/davidemms/OrthoFinder
Question: Dear lazyweb, here is the new menu for `samtools index`:

```
Usage: samtools index [-bc] [-m INT] [out.index]
Options:
-b Generate BAI-format index for BAM files [default]
-c Generate CSI-format index for BAM files
-m INT Set minimum interval size for CSI indices to 2^INT [14]
```

What is the **CSI** index compared to **BAI**? Should I better use it instead of the 'old' BAI index? What's the benefit?

Answer: CSI also allows to choose a compromise between the size of index and region query performance (the `-m` parameter). The less the interval size, the larger the index is, but for region queries less computational work is required, because the index is more fine-grained.

BAI is a special case of CSI, with some minor differences.
Question: Hi, I want to use bedops to analyze some of the .bed files from ChIP-seq, but one of the files I can't go through the analysis, even can't sort-bed. It keeps saying the BED row length exceeds capacity at line 1. It doesn't help even after I deleted the first row of my .bed file.

```
unknownb8f6b1106ced:chip-seq fuxiaoyong$ sort-bed era_inpegf_mcf7_n3_hg18_f20_nr.bed
BED row length exceeds capacity at line 1 in era_inpegf_mcf7_n3_hg18_f20_nr.bed.
Check that you have unix newlines (cat -A) or increase TOKENS_MAX_LENGTH in BEDOPS.Constants.hpp and recompile BEDOPS.
```

I am very new for using tools to analyze NGS data, even for the Mac OS X system. So, please help me and your detailed explanation will be greatly appreciated!!

Answer: The issue is probably not blank lines, which the current version of `sort-bed` removes, but a lack of correct newline characters, since Excel (even the Mac version) adds Windows-specific line delimiters to the text files it exports.

If the file doesn't end with a UNIX newline character, then you will get an error message. The usual way to fix this is to use `dos2unix` on the output from Excel, or to avoid handling data in Excel where possible.

If clearing out empty lines fixes sorting for you, then look at upgrading your version of BEDOPS and `sort-bed` to the current release (presently v2.4.2), as you might be running an older binary.
Question: nhmmer is the DNA pipeline version of hmmsearch within the HMMER software. Recently, I've been comparing two paired-end Illumina libraries with the following:

1. Trimmomatic
2. bwa aln, forward/reverse trimmed reads to database of known reference genes
3. Parsing the SAM file for mapped reads (samtools -F 4) and analyzing those

I ran the same **raw read files** through nhmmer and analyzed the intersection of the alignment and HMMER pipelines. When I look at reads that aligned but were not called by HMMER, bit 16 (reverse strand) is set on virtually all of them. Additionally, when the read information is pulled from the SAM file (where it has been appropriately reverse complemented) and re-run through HMMER, nhmmer does find ~70% of them.

So it appears that HMMER doesn't consider the reverse complement of the database being queried, though it does consider the reverse orientation when the query is of the same strandedness with respect to the database. This is in contrast to the manual, which claims that it does search both strands of the database, but this doesn't explain the findings here.

Has anyone had similar experiences, or does someone have inside knowledge on how or if nhmmer considers reverse complements?

Answer: HMMER does in fact do this correctly; it was an issue with an intermediate step written by myself. Lesson learned; do more internal testing on code.
Question: Hello everyone.

I need a code that can give me chunks of secondary structure elements from a given RNA structure and its sequence. e.g. if I have a sequence and its structure like this:

sequence: UGCCCUAGUACGAGAGGACCGGAGUG

Structure: `...(((........))).........`

I should be able to get the output like:

```
Loop1: UGC
Loop2: AGUACGAG
Loop3:GGAGUG
Stem1: CCUAGG
```

The position of every dot corresponds to the position of it respective nucleotide. I have written this code but it gives me all the nucleotides which form loop as one string and all those forming stems as one string. I want fragments.

Answer: [forgi][1] ([github page][2]) is a Python package for working with RNA secondary structures. I've used it in the past for tasks like this.

Here's an example. By the way, might want to check your expected output against the example sequence you provided.

```
import forgi.graph.bulge_graph as fgb
bg = fgb.BulgeGraph()
bg.from_fasta(
"""\
>demo
UGCCCUAGUACGAGAGGACCGGAGUG
...(((........))).........
""")
for d in bg.defines:
print d, bg.get_define_seq_str(d)
```

Here's the output. Notice it keeps the two sides of the stem separate, and the 5' and 3' "stems" are listed separately.

```
f1 ['UGC']
h0 ['AGUACGAG']
s0 ['CCU', 'AGG']
t1 ['ACCGGAGUG']
```

Edit: [forna][3] by the same author is useful for visualizing structures and is helpful for debugging and troubleshooting this sort of analysis.

Edit 2: In response to the comment asking about providing a string variable to `from_fasta`:

It's not documented, but after playing around it looks like this method requires 1) the record must start with ">" -- not even a newline is allowed, and 2) a trailing newline at the end of the structure line. So if you have your strings `sequence` and `structure`, you can build a string that `from_fasta` will accept like this:

```
text = '\n'.join(['> ', sequence, structure, '\n'])
bg.from_fasta(text)
```

or if you want to include a sequence ID:

```
text = '\n'.join(['\n', '> ' + seqid, sequence, structure, '\n'])
bg.from_fasta(text)
```

[1]: http://www.tbi.univie.ac.at/~pkerp/forgi/index.html
[2]: https://github.com/pkerpedjiev/forgi
[3]: http://nibiru.tbi.univie.ac.at/forna/
Question: I need to check which of the predicted de novo transcripts have polyA at 3'. Are polyA tails filtered at some step of tophat/cufflinks processing or this information is retained? If this info is retained, where do I find it. thank you

Answer: Most reads from bridge amplification wouldn't contain polyA's nor would they ever make it past alignment. I think you are confusing de novo transcriptome assembly with ab initio isoform discovery.
Question: Can anyone explain what the allowed values for the SAM header (@HD) VN tag are, and what they imply? For example, the [current specification PDF][1] shows:

@HD VN:1.5 SO:coordinate

What is implied by stating your SAM file is versioned at "1.5", and do any tools actually care about this value? I'm asking because I was about to create an unaligned SAM file from scratch, and am at a loss about how to responsibly set this value.

[1]: http://samtools.github.io/hts-specs/SAMv1.pdf

Answer: Ah the joys of undocumented things. For what it's worth, most tools completely ignore the version number (samtools sort will even just use VN:1.3 if none is specified in the file). I think the only one to watch out for is 1.4 and later, since that added CIGAR operations...though almost no one uses them.