cannin/biostars_qa huggingface dataset:
Question: In SNPEFF-IMPACT tool, What is difference between STOP_GAINED and NON-SYNONYMOUS_STOP?

In SNPEFF, effects are categorized by their impact including (High, Moderate, Low and Modifier).
I would like to know, what is difference between STOP-GAINED and NON-SYNONYMOUS-STOP?

STOP-GAINED has been categorized in High impact and NON-SYNONYMOUS-STOP has been categorized in Low impact.

Reference : SNPEFF table

Answer: Without looking at it, my guess is that the variant changes the original stop codon to a different stop codon (that is, it switches among TAG, TAA, and TGA). So the codon changes (it's non-synonymous), but the effect for the gene should not.

Edit: What I described is a SYNONYMOUS_STOP. As I understand it, a NON_SYNONYMOUS_STOP can't arise, and exists in the code only as an oversight. As evidence of this, only a SYNONYMOUS STOP is explicitly documented in the FAQ. Furthermore, in one of my datasets, I observe over 1000 SYNONYMOUS_STOP calls and no NON_SYNONYMOUS_STOP calls.

More explanation from the relevant code:

if (aaOld.equals(aaNew)) {
// Same AA: Synonymous coding
if ((codonNum == 0) && codonTable.isStartFirst(codonsOld)) {
// details removed
} else if (codonTable.isStop(codonsOld)) {
if (codonTable.isStop(codonsNew)) effectType = EffectType.SYNONYMOUS_STOP;
else effectType = EffectType.STOP_LOST;
} else effectType = EffectType.SYNONYMOUS_CODING;
} else {
// Different AA: Non-synonymous coding
if ((codonNum == 0) && codonTable.isStartFirst(codonsOld)) {
// details removed
} else if (codonTable.isStop(codonsOld)) {
if (codonTable.isStop(codonsNew)) effectType = EffectType.NON_SYNONYMOUS_STOP;
else effectType = EffectType.STOP_LOST;
} else if (codonTable.isStop(codonsNew)) effectType = EffectType.STOP_GAINED;
else effectType = EffectType.NON_SYNONYMOUS_CODING;
}

There is no way for the old and new codons to be different (to get into the outer else block) and for both codons to be stop. If they're both stop codons, they'll be the same and go to the first if branch. So the NON_SYNONYMOUS_STOP branch is never reached.
Question: Hi,

I have a number of atacseq bams and their total read numbers differ. I did peak calling with and without --SPMR for normalization and there is no difference between the two ......broad_peaks.broadPeak files. yet there is a difference in the resulting two ....broad_treat_pileup.bdg files. While using callpeak, does --SPMR have no effect on the signalValues? If that is the case, should I do another peak calling from the bdg files?

My end goal is to see what happened to some specific regions, whether they got condensed or decondensed and I cannot make sure whether this read depth difference effect the signalValues.

Answer: Whether you use the SPMR flag or not, MACS2 will normalise the data internally to call peaks. The SPMR flag only effects the signal track produced. With the flag present the signal will be normalised to reads per million, this is for comparison with other samples which have been sequenced to different read depths.
Question: Here is an example of RNA velocity from [Manno .et .al 2018](https://www.nature.com/articles/s41586-018-0414-6).

![enter image description here][1]

Colour of cells denote cell type. Layout denotes how cells cluster by some measure of similarity. The direction of the arrows denote the future state of the current cells. Or I guess one can say the differentiation trajectory.

What do length of the arrows denote? Can we say something about the propensity of that cell to differentiate in that direction? Can we say if there are more or less spliced or unspliced transcripts?

My understanding is that an arrow is based on the number of spliced and unspliced transcripts for each gene, summarised over all genes and further summarised over all cells in a small neighbourhood.

Other comments on understanding RNA velocity plots are welcome. I am trying to understand this better.

[1]: https://i.ibb.co/YDmQJFN/41586-2018-414-Fig4-HTML.png

Answer: This paragraph in Svensson and Pachter, *Molecular Cell* 2018 (https://www.sciencedirect.com/science/article/pii/S1097276518307974) might help understand the concept of these arrows:

> "Long arrows, corresponding to large changes in gene expression, will
> in the low-dimensional representation imply that cells are undergoing
> rapid differentiation, while short arrows can be intuitively
> associated with terminally differentiated cells maintaining
> homeostasis. While some care needs to be taken when interpreting the
> lengths of arrows in low-dimensional representations, RNA velocity
> representations can be biologically informative. In principal
> component analysis, the vector lengths are meaningful, scaling
> linearly with RNA velocity estimates in high dimensions. We note that
> this will not be the case in t-SNE representations (Van der Maaten and
> Hinton, 2008), though directions remain meaningful and can be
> interpreted as “RNA arguments” in vector nomenclature."
Question: what is the biological meaning of chromosome cytoband except centromere? they contain more genes or more active than other regions?

Answer: The cytobands/karyotype bands/chromosome bands/G banding etc are old fashioned data. They date back to before we had genome sequences, when the way we used to identify chromosomes and chromosome regions by staining metaphase chromosomes with Giemsa, then viewing them with a microscope. Heterochromatic (closed chromatin, low gene density) regions stain darker than euchromatic (open chromatin, high gene density) regions. Bioinformatically, they're obsolete as genome coordinates are far more accurate and meaningful, however people still use them as a shorthand for identifying genomic regions; it's much quicker, easier and more memorable to say 14q21.3 than it is to say 14:46695396-50395063. I wouldn't use them for looking at heterochromatin vs euchromatin these days either, I would instead look at actual gene density and publicly available DNase sensitivity for open/closed chromatin.
Question: Hello community!

I just started with bioinformatics, so I hope you can help me with understanding some basic things. I have RNA-Seq 100bp pair-ended reads from the *Lepidoptera *insect, for which there is no genome available. In total there were 2 treatments and one control, so for the assembly I merged all of them into two files (R1 and R2). So the goal is to make an assembly, and then continue with mapping, counts and statistic.

For the assembly I used next tools: Trinity, Velvet/Oases, Bridger and CLC assembler. From the last one I just got the ready-to-go assembly, to which I compared other assemblers.

1. Trinity was executed with default parameters (k=25),
2. Velvet: velveth/velvetg k=21,51,2; For each assembly run oases and compared n50 and total number of contigs. Chose the range of contigs from 25 till 39, with kmer=35, velveth/velvetg on the oases results and then merged them with oases.
3. Bridger. Default parameters with kmers 25 and 27.

After I compared assemblers with Quast. As a result got the next table

```
Oases Trinity Bridger 25 Bridger 27 CLC
#contigs 10706 15579 14881 14632 19985
#contigs (>= 0 bp) 36576 40429 37011 36162 35589
Largest contig 27732 18062 27729 27887 40347
Total length 12451630 19701362 20030913 19559912 36090599
N50 1324 1528 1653 1651 2584
# N's per 100 kbp 0 0 0.05 0.21 212.61
# N's 0 0 11 41 76733
```

So, I could not reach the N50 of CLC, by using default parameters. But at least I can see how do they reach it. The number of mismatches in CLC is super high, so I was wondering - is it really worth to increase the total length and n50 by putting double amount of N's into the contigs?

And the next question - for some reason I though it would be a good idea to increase the input sequence length by merging R1 and R1, R2 and R2 reads. So basically doubling the amount of data assembler should deal with. So far I tested only Bridger (k=25), but the results are surprising, N50 increased to 1720, and the total length and the largest contigs became almost the same as in CLC results. What I think is happened that during the building of the graph low quality k-mers had a better chance to survive, thus increasing the contig length.

Thanks in advance for any kind of responds!

Answer: Assembly is somewhat of a "dark art" and that, as I found, means trying out many different methods/parameters until one seems to work. ... not the most scientific approach though. Merging the data if it is from different measurements and conditions should be just fine, even recommended. It is possible that some transcripts get only expressed under one of the conditions.

Fundamentally the problem is that N50 is not a good measure of assembly quality (although it usually correlates with it) and is the only one you can easily compute right away.

What I would recommend is to align assemblies against each other. You may find good agreement for some transcripts, these are more likely to be correct, then many disagreements and you can perhaps find some patterns.
Question: I have a file with 2 columns, the first column have a read ID from fastq, the second column have a ncbi taxi ID, like this:

ST-E00243:601:HWJFHCCXY:8:1101:19431:10468	9606
ST-E00243:601:HWJFHCCXY:8:1101:19451:10468	14230
ST-E00243:601:HWJFHCCXY:8:1101:19471:10468	10468
ST-E00243:601:HWJFHCCXY:8:1101:19492:10468	1512
ST-E00243:601:HWJFHCCXY:8:1101:19512:10468	512
ST-E00243:601:HWJFHCCXY:8:1101:19532:10468	1421067

and I have an archive with the taxid list of a specific taxos, like this:

2485233
2485231
2059665
2029516
2022430
2022429
1987726
1980986
1738445
1737346

I want to extract the lines with the taxid present on archive 2, but, some of the taxids presents the same number present on read ID. So I tried extract only the correct lines with fgrep, but, fgrep can't specify a column (in the case the column 2).

Anyone can help me? Briefly I need extract the whole line detecting a pattern in column 2 using a file with patterns.

Thanks !

Answer: Assuming tab-delimited files:

cat file1.txt
2485233
2485231
2059665
2029516
512
2022430
2022429
9606
1987726
1980986
1738445
1737346

cat file2.txt
ST-E00243:601:HWJFHCCXY:8:1101:19431:10468	9606
ST-E00243:601:HWJFHCCXY:8:1101:19451:10468	14230
ST-E00243:601:HWJFHCCXY:8:1101:19471:10468	10468
ST-E00243:601:HWJFHCCXY:8:1101:19492:10468	1512
ST-E00243:601:HWJFHCCXY:8:1101:19512:10468	512
ST-E00243:601:HWJFHCCXY:8:1101:19532:10468	1421067

awk 'FNR==NR {a[$1]==$1; next} {if ($2 in a) print $0}' file1.txt file2.txt
ST-E00243:601:HWJFHCCXY:8:1101:19431:10468	9606
ST-E00243:601:HWJFHCCXY:8:1101:19512:10468	512
Question: Hi guys,
I have a question about how to choose a reasonable threshold to define genes to be expressed (not differentially expressed but simply expressed) to compare gene expression measurements from different platforms. In other words I have to check if my list of genes is expressed using an RNA seq dataset, using a data set of one-color microarrays and another of two colours microarrays. If for RNA Seq is generally "accepted" (but I know that it depends on the experimental condition) an rpkm >1 I don't know how to set a threshold for microarrays. Any suggestion about this?

Thank you in advance

Answer: There is no generally accepted threshold in RNA-seq. You can come up with a threshold using zFPKMs, but even that is a bit hand-wavy. For microarrays, you generally check either that the PM signal is sufficiently above the MM signal and/or that a probe is sufficiently above that of the background control probes. Not that these are sort of apples and oranges measurements, in that they're really not comparable.
Question: I am reading two papers of Bolstad about quantile normalization.

In "Probe Level Quantile Normalization of High Density Oligonucleotide Array Data" he writes about quantile normalization:

> One problem with this method is that in the tails in particular,where we might expect greater differentiation between chips,the normalized values are going to be identical

However, in "A comparison of normalization nethods for high density oligonucleotide array data based on variance and bias" he writes:

> This would be most problematic in the tails where it is possible that a probe could have the same value across all the arrays.

If I understood it correctly, in the first one he states that in tails (say where we have high intensities) the probe will have greater differentiation between the chips, so the values of the probe will not be the same between the chips.

In the second he states that a probe still could have the same value across the chips.

Isn`t it controversially? Or I understood the first paper wrongly and he is talking about the probeset in the first paper?

edit: or in the first paper he talks about unnormalized values and in the second about already normalized?

Thanks in advance.
Tonja

Answer: He is referring to the same effect of quantile normalization in both sentences. When he says "greater differentiation" he means in reality, before quantile normalization. When he says "same values across all arrays" he is referring to post-normalization values. Consider an extreme case of a probe intensity for an interesting probe that moves to the 100th percentile in a small subset of samples representing a treatment group. After normalization, that probe's intensity in that subset of samples will be identical to the 100th percentile intensities of all the other arrays, potentially dampening an interesting change. Or, if a single probe (corresponding to a very abundant mRNA) was at the 100th percentile of all samples, after normalization it will have a single, identical, constant value across the entire dataset.
Question: What are the characteristic features of Chip-seq and RNA-seq data? If I have fastq files which are the results of Chip-seq and RNA-seq experiment is it possible to differentiate then, for example by comparing to Chip-seq input result, which is explicitly marked?

Answer: **RNA-seq** : (In eukaryots) splicing (some reads must be split to map), very uneven read coverage, especially in total RNA-seq where rRNA reads dominate.

**ChIP-seq** : relatively even coverage in the input fraction.

To check genome coverage requires mapping, but over-represented sequences can be analysed with [fastqc][1] and blast to provide a quick indication directly from the fastq files : If there are over-represented sequences that corresponds to highly expressed genes, then you are dealing with RNA-seq data.

[1]: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Question: Hi all,

I would like to filter / annotate my VCF and/or BAM files using dbSNP as part of the input during the steps of the variant calling workflow (other than the base recalibration step), are there any steps along the workflow where one can do so?

Or are you using Annovar or a similar program for that once you have a filtered VCF file?

Our workflow is based on the following: http://www.htslib.org/workflow/#mapping_to_variant

In addition to the variant caller mentioned there we also use FreeBayes & GATK for variant calling.

Thank you!

Answer: HI,

You can use [SnpEff][1] to annotate and [SnpSift][2] to filter variants using dbSNP.

Some useful information [here][3] (example 5).

[1]: http://snpeff.sourceforge.net/
[2]: http://snpeff.sourceforge.net/SnpSift.html
[3]: http://snpeff.sourceforge.net/protocol.html#exFilter
Question: Very basic question here (so that I'm afraid to ask). Suppose I have 50 pairs of data and I want to linearly model my data and predict a single new pair to my single input. What should I do?

a= matrix (rnorm(100),50)
l= lm (a[,2]~a[,1])
predict (l, 2.1)

doesn't work.

Answer: Try this example:

# reproducible example input
set.seed(1); a <- data.frame(matrix(rnorm(100), ncol = 2))
head(a)
# X1 X2
# 1 -0.6264538 0.3981059
# 2 0.1836433 -0.6120264
# 3 -0.8356286 0.3411197
# 4 1.5952808 -1.1293631
# 5 0.3295078 1.4330237
# 6 -0.8204684 1.9803999

# get fit
l <- lm (X2 ~ X1, data = a)

# new input data with only X1 column, because we want to predict X2
newData <- data.frame(X1 = 2.1)

newData$X2 <- predict(l, newData)
newData
# X1 X2
# 1 2.1 0.02624987
Question: Physical (or fragment) coverage means counting also the unsequenced part of the fragment between properly mapped read pairs. As output I'd like to get a frequency table that tells me how many bases are covered how many times. I could not find an easy solution for this problem using standard tools only (SAMtools, BEDtools, GATK).

EDIT: The reason why I am interested in physical coverage is because I think this is quite a useful metric when it comes to assess the power of a WGS data set to detect structural variations.

Answer: You can do this in three steps using bedtools. First convert your paired-end BAM file to a BEDPE formated file of fragments with bedtools bamToBed, then use cut to convert BEDPE to BED, then use bedtools genomecov to compute coverage from the resulting BED file of fragments. See this the section after NOTES in this post at the bedtools mailing list for more information: https://groups.google.com/forum/#!msg/bedtools-discuss/kulRA_7Ybow/4Be6otCybZgJ
Question: Hi

I need to design primers for around 40,000 sequences. After doing this task with Primer3, I found that it took a very long time.

I I tried to accelerate primer3 operation with GNU parallel, but I cannot managed to successfully use GNU parallel to split input file and do multi-thread operation. Somehow primer3 still ran on 1 core only.

My command is as the following:

cat fasta.p3in | parallel --round-robin -j 12 --pipe --recend "=" /Tools/primer3/primer3-2.3.6/src/primer3_core > fasta.p3out

Could anyone tell the correct way to use GNU Parallel along with Primer3? Thanks a lot!

Answer: Your biggest mistake was probably that your records contain `=` on every line, but only `\n=\n` is a record separator. Using the command `wc` or `--files cat` is great for debugging that kind of problems.

Your second mistake is that `--block-size` defaults to 1M: So the first instance may simply gobble up everything.

This ought to work (untested, as I have neither access to `fasta.p3in` nor to `primer3`):

cat fasta.p3in | parallel -N1 --round-robin --pipe --recend "\n=\n" --cat /Tools/primer3/primer3-2.3.6/src/primer3_core > fasta.p3out

You can possibly leave out `--cat` if primer3 reads from STDIN. If GNU Parallel takes up significant time, increase `-N1`: With 40000 records it is probably OK to split on bigger chunks than 1 record.
Question: Hi,

I have an expression profile dataset with a list of genes and normalized expression data. What i want is to only look at genes on the Y chromosome.

Most of the gene annotation packages with R or otherwise look for gene ids/probe-set ID to map to chromosomal information. Can you please help me out with a way to extract chromosomal information ( preferably through an R package) with only the gene names ?

Thanks !

Answer: The trick is to use the "org" package for your organism. Here I show human, but you can substitute....

library(org.Hs.eg.db)
symbols = c('BRCA1','TP53','GAPDH')
select(org.Hs.eg.db,keys=symbols,cols=c('SYMBOL','CHR'),keytype='SYMBOL')

Which will return:

SYMBOL CHR
1 BRCA1 17
2 TP53 17
3 GAPDH 12
Question: Simple question.... Has anyone produced a plugin for adding [ExAC][1] allele frequencies to Variant Effect Predictor (VEP) annotation?

[1]: http://exac.broadinstitute.org/

Answer: Available for a while, but I forgot to post it: https://github.com/ensembl-variation/VEP_plugins/blob/master/ExAC.pm
Question: Hi all,

does anyone know how to check if a tree has been correctly converted to ultrametric with ete3 in python? Kinda like the is.ultrametric() function of the R ape package?

Thanks!

Answer: try something like this:

from ete3 import Tree

def is_ultrametric(tree):
last_distance = None
distance = 0
for post, node in tree.iter_prepostorder():
if post:
distance -= node.dist
else:
if not node.is_leaf():
distance += node.dist
else:
if last_distance is None:
print(last_distance)
last_distance = distance + node.dist
elif last_distance != distance + node.dist:
return False, last_distance
return True, last_distance

t = Tree('((a:1,(b:0.5, d:0.5):0.5):1,c:2);')
print(is_ultrametric(t))

t = Tree()
t.populate(10, random_branches=True)
print(is_ultrametric(t))
Question: I have several assembled, annotated transcriptomes of non-model organisms (fish). I would like to identify what transposons are present in these transcriptomes (LTR, SINEs, LINEs etc) and compare them. I am currently trying to use RepeatMasker to identify them using a clade specific library from RepBase but I was thinking of trying to look for them de novo with REPET. Just looking for suggestions and ideas really. Thank you.

Answer: Tandem Repeat Finder is also a good option for finding denovo repeats.

You can also try, Censor for finding different kinds of repeats from your sequences. Check this paper analysing transposons using censor and repeatmasker.
Question: Hi,

I have a pipeline where I would previously go on ensembl biomart and retrieve data but would now like to automatize a bit more and get my data from biomart directly in a Python script using the Bioservices package. My initial input is a list of RefSeq mRNA identifiers and then I use biomart to retrieve the corresponding ensembl gene identifiers. This works fine. Then I want to give these ensembl identifiers to biomart and to get the corresponding CDS sequences, accompanied by the ensembl gene ID, the ensembl transcript ID and the chromosomal start positions of exons*. At this point, I expect to get back a bit less than 2000 sequences (because this is what I get when I do this manually through my browser). In actuality, I only get 8 sequences and they are accompanied not by ensembl IDs but rather by LRG IDs! The sequences I do get look about right though. If anybody had any idea what was going on here, I would really appreciate your help. I am using Python3.4 on Mac OS X 10.9.5. Here is my code:

```
input_file_name = "single_exon_pcgenes_refseq2.txt"

from bioservices import *
import re
e = ensembl.Ensembl()

bm = BioMart()
bm.datasets("ensembl")
bm.add_dataset_to_xml("hsapiens_gene_ensembl")
refseq_ids_file = open(input_file_name)
refseq_ids = refseq_ids_file.readlines()
refseq_ids_file.close()
refseq_ids = [re.sub("\n","",i) for I in refseq_ids]
refseq_ids = ",".join(refseq_ids)
bm.add_filter_to_xml("refseq_mrna", refseq_ids)
bm.add_attribute_to_xml("ensembl_gene_id")
xml_query = bm.get_xml()
ens_ids_raw = bm.query(xml_query)

bm.new_query()
bm.add_dataset_to_xml("hsapiens_gene_ensembl")
bm.add_filter_to_xml("ens_hs_gene", ens_ids_raw)
bm.add_attribute_to_xml("ensembl_gene_id")
bm.add_attribute_to_xml("ensembl_transcript_id")
bm.add_attribute_to_xml("exon_chrom_start")
bm.add_attribute_to_xml("coding")
xml_query_seq = bm.get_xml()
raw_sequences = bm.query(xml_query_seq)
```

This is one of the four sequences contained in `raw_sequences` at the moment (I deleted the middle part of the sequence to make it easier to read and added newlines):

```
ATGGAGGGGATCAGTATATACACTTCAGATAACTACACCGAGGAAATGGGCTCAGGGGACTATGACTCCAT
AACCTCTACAGCAGTGTCCTCATCCTGGCCTTCATCAGTCTGGACCGCTACCTGGCCATCGTCCACGCCA
GAGTCTTCAAGTTTTCACTCCAGCTAA LRG_51 LRG_51t1 5001;7244
```

*This might seem like a strange way of doing things but I can't ask for the sequences directly in the first step using the RefSeq IDs as input because this would only give me the sequences of the particular transcripts whereas I need the sequences of all the transcripts from each corresponding gene. But I need to filter them first to make sure that each of those genes has at least one single-exon transcript, which is why I need to go through the transcript IDs in the first place (I get those from UCSC).

Answer: I got it!

Sorry, I should have kept on playing with it for a few more hours before turning to you guys but it was starting to look a bit grim there a while ago.

It turns out I was simply using the wrong filter name. So

bm.add_filter_to_xml("ens_hs_gene", ens_ids_raw)

should have been

bm.add_filter_to_xml("ensembl_gene_id, ens_ids_raw)

Sorry, that was a dumb one!