cannin/biostars_qa huggingface dataset:
Question: I have noticed a long time ago that the true read length may not exactly match the promised read length in the analysis name e.g. 2x100bp Illumina data frequently consists of 101bp reads. Today I looked into a dataset that systematically has 97 or 96 bp reads for the promised 2x100bp across many files. Why does this happen and should it be interpreted in any way?

Thanks in advance for the comments!

Answer: The explanation for having one base more than expected (e.g. 101 instead of 100), is that for Illumina data, the phasing for each position in the read depends also on the data from the next base, so the last base in the read will always have a much worse base quality, because there is no data from a next base. So by having 101 instead of 100 bases you should get 100 bases with reliable base qualities.

For the cases where there are 96 or 97 bases, I'm guessing that either the run ran out of reagents, or else the base qualities at the 3' ends of the reads were poor and someone has already trimmed the last few bases.
Question: Hi All,
I am working on Single-cell data and I am using Seurat for the data analysis. I want a Violin plot showing relative expression of select differentially expressed genes (columns) for each cluster as shown in the figure (rows) (all Padj < 0.05).

![Violin Plot][1]

It would be really helpful if you can let me know how to plot as shown in the image.

Thanks a lot for your kind help,

Dave.

[1]: https://i.ibb.co/YkR3jSF/Screen-Shot-2020-07-27-at-8-27-44-PM.png

Answer: Hi,

The above plot isn’t from Seurat AFAIK, it is from `scanpy()`.
Last time I checked this multi gene plot function wasn’t available but one can play around and generate it. Ming has a blog post with code around `ggplot()`. If you are not familiar with `python` and underlying packages of single cell there then with `Rstats` and `Seurat()` along with other packages you can still produce it. Please check the below link to have the code snippet and try it out at your end. Good luck.

[Multi gene violin plot][1]

[1]: https://divingintogeneticsandgenomics.rbind.io/post/stacked-violin-plot-for-visualizing-single-cell-data-in-seurat/
Question: I have 2 sets of bed files from ChIP-seq, A and B, and a certain number of peaks overlap. Now I want to put some kind of enrichment score / p-value on that overlap - if only 2 peaks overlap and A and B contain both 10 000 peaks, it's not really significant, while if 90 peaks overlap, A is 10000 peaks and B is 100, it's more significant.

The classical way to do this would be to do a hypergeometric test - all the variables of the test are easily filled in (number of draws, number of successes, number of 'failures') except one, the total population number. In this scenario it would be something similar to 'the total number of peaks one could potentially draw from the genome', which is impossible to estimate, and could be biased (regions that don't chip well, unmappable regions, etc).

What's the best way to put a 'score' on an overlap between A and B then? The ultimate goal is to display some kind of heatmap, as I'm doing the comparison between many A's and many B's.

Answer: Something to this effect is available in the latest release of bedtools as [bedtools fisher][1].

It uses bases, rather than counts of peaks. Currently, it doesn't allow you to exclude regions like centromeres or "quiescient" regions that are unlikely to harbor overlaps and I have a PR with that as a work in progress, but it won't be ready for a while.

I have a ticket open to allow the use of counts of intervals, but it will suffer from determining the denominator as you suggest. I think that the total number of possible intervals could be the number of times an interval of the average size could fit into the genome.

For your stuff, you may want to check out [this paper][2] and the associated software.

[1]: http://bedtools.readthedocs.org/en/latest/content/tools/fisher.html
[2]: http://www.ncbi.nlm.nih.gov/pubmed/22262674
Question: Hi, I have some doubts on cpm function in edgeR package. Given this initial part of the code:

> library(edgeR)
> data_clean<- read.table("unito.txt", header = TRUE)
> group<-factor(c("A", "A", B","B"))
> y <- DGEList(counts=data_clean,group=group)
> y <- calcNormFactors(y)

I obtain different values if I calculate cpm in this way

> cpm <- cpm(y, log = FALSE, normalized.lib.sizes=TRUE)

Or in this one

> cpm <- cpm(y$counts, log = FALSE, normalized.lib.sizes=TRUE)

Why are them different? What does it change? What is the correct way to calculate them? I noticed that if I calculate them with the standard formula "by hand", I obtain the same values as the second formula (y$counts). In the edge manual it seems just to use the y argument.

Thank you,
Francesca

Answer: This is normal and expected. If you use the first chunk the function will use the `norm.factors` from `calcNormFactors` to correct the counts for both library size **and** composition, see [this][1] video for the theory of library normalization.

With the second chunk of code you are providing only the count matrix so the function has no `norm.factors` therefore only correct per-million so for library size. That is why this matches with your `by hand` calculation.

The first method is the one you should use.

[1]: https://www.youtube.com/watch?v=Wdt6jdi-NQo
Question: I am doing a school assignment and I find a snippet of code extremely confusing. I have copied and pasted it from the assignment to here.
library(limma)
library(tidyverse)
library(ggplot2)
library(knitr)
library(Biobase)
library(GEOquery)

eset <- getGEO("GSE4051", getGPL = FALSE)[[1]]

This part is what I find confusing. Later on pData(eset) gets called and the value of the method being called has been changed. What is going on? pData(eset) is actually an object or something? I don't understand why this isn't being saved as a global variable.

pData(eset) <- pData(eset) %>%
mutate(sample_id = geo_accession) %>%
mutate(dev_stage = case_when(
grepl("E16", title) ~ "E16",
grepl("P2", title) ~ "P2",
grepl("P6", title) ~ "P6",
grepl("P10", title) ~ "P10",
grepl("4 weeks", title) ~ "4_weeks"
)) %>%
mutate(genotype = case_when(
grepl("Nrl-ko", title) ~ "NrlKO",
grepl("wt", title) ~ "WT"
))

pData(eset) %>%
with(table(dev_stage, genotype))

Answer: `getGEO(...)` returns an object, and you are assigning this object to the `eset` variable. You can check the object class using `class(eset)`. The object type is likely S4, which for us means it has various slots that store different pieces of information.

Whatever object type eset is likely had "setter" and "getter" generic functions defined for a particular slot. The getter `pData(eset)` extracts data from an object slot, and the setter `pData(eset) <-` adds data to an object slot.

`pData(eset) <- pData(eset) %>% ...` is getting data from an object slot, manipulating the data, and then saving it back into the same object slot.

You can read more information about this topic in the [S4 objects chapter of Advanced R](https://adv-r.hadley.nz/s4.html).
Question: **Dear Biologists and Bio-IT-specialists,**

based on your knoweledge and experience I'd like to know how many variants should be discovered from .BAM file (final result, after annotation and filtration)?
This question is very general, so I'll specify it a bit.

Let's take some NGS data (fastq) from ENA databse. I performed whole analysis and here's my results:

**WES data, paired end, Illumina HiSeq 2000, material - cancer cell line. Software that I used:**

1. FastQC - (before and after trimming adapters)
2. Trimmomatic - trimming adapters
3. Subread - mapping to genome - hg19
4. Bamtools - filter, sort, index
5. FreeBayes - variant discovery
6. SnpSift i snpEff - filter and annotation of discovered variants.

After annotation of variants I run:

(QUAL > 1) & (QUAL / AO > 10) & (SAF > 0) & (SAR > 0) & (RPR > 1) & (RPL > 1)

as a filtering arugemnt of SnpSift.
As a result I have .vcf file with **32866 variants**. Do you think that this number might be correct for mentioned above criteria?

Why do I ask such question like this?

A few days ago I got several .vcf (ready-to-use) files of WES experiment of another cancer but from patients... Analysis were performed by some specialist that I don't know. I was a bit confused becaue each .vcf file contains **300000-350000** variants with filter **PASS**. To be honest I don't know the pipeline - of analysis that returned above results, however variant discovery was performed by SAMTools...

**So my questions are:**

Do you think that my analysis and tools that I used are correct? I'd like to avoid false negatives when some mutation is present and it's absent in my result .vcf file...

How many variants of WES experimrnt (final result) analysis do you usually have?

Best regards,

Adam

Answer: The total number of variants really depends on the sample. For cell lines, as they are often in culture for decades, a plethora of mutations might have been acquired over the years which do not really have something to do with the underlying disease. IMHO, the only use of cell line variants is to check if they are suitable models for experiments in case that mutations from patients are also present in the line. Variant number also depends on the disease, as there are cancer types such as AML with low mutation rates, while others such as lymphoma and most skin- and lung cancers have high mutation rates. That having said, whatever the numbers in your cell line are, I think it is not representative for anything. Also keep in mind that most of these mutations are probably germline variants, which you cannot discriminate from somatic ones due to the lack of a matched-normal control. What you can do is to filter out common variants in the human population, e.g. removing all variants that are reported by the 1000Genomes project. The Ensembl Variant Effector Predictor (VEP) can do that for you. But as I said, cell line mutations must be considered carefully and cannot replace primary material.
Question: Hello,

I'm wondering whether there is a program that could calculate chromosome sizes from any fasta file? The idea is to generate a tab file like the one expected in bedtools genomecov for example.

I know there's the fetchChromSize program from UCSC, but not all genomes are available over there (I need TAIR10 for instance). I've read this topic already.

I would like a tool that can deal with any genome regardless of the database. If it doesn't exist I guess it's possible to just parse fasta files, but I'd be surprised if no one else had done it before!

Cheers

Answer: pip install pyfaidx
faidx input.fasta -i chromsizes > sizes.genome

That's what you're looking for if you want to use `bedtools genomecov`, but you can transform to a BED file as well using `-i bed`.
Question: When I search a gene("METRNL") in UCSC, I found the gene located in "chr17" and "chr17_KI270860v1_alt" (see: goo.gl/xzz2NA). for hg38.

What is the relation between these two?

Answer: Copying from this announcement about hg20/GrCh38:

- **Alternate sequences** - Several human chromosomal regions exhibit sufficient variability to prevent adequate representation by a single sequence. To address this, the GRCh38 assembly provides alternate sequence for selected variant regions through the inclusion of alternate loci scaffolds (or alt loci). Alt loci are separate accessioned sequences that are aligned to reference chromosomes. This assembly contains 261 alt loci, many of which are associated with the LRC/KIR area of chr19 and the MHC region on chr6.

The same principle applies to some regions in hg18/19, etc.
Question: Hi All,

I am using blastn to blast my data against tRNAs which is extracted from RFam database. Then I picked one match from the blast result to confirm it by directly searching my read in RFam online (http://rfam.sanger.ac.uk/search).

Here is my read which have a 100% match with a tRNA:
CCCATTCTTGCGACCCGGGTTCGATTCCCGGGCGGCGC

This is the matched tRNA from RFam:
GCGCCGCTGGTGTAGTGGTATCATGCAAGATTCCCATTCTTGCGACCCGGGTTCGATTCCCGGGCGGCGCA

Clearly it's a perfect match (in bold). But when you search my query in RFam, there is no hit. However, the subject tRNA can be found in RFam.

Could anybody tell me what's the problem?

Answer: I think that it is because your tRNA sequence is truncated by having short reads. Covariance models use (among other things) information about the internal base-pairing of the RNA, in this case the cloverleaf structure.

(By Yikrazuul (Own work) CC-BY-SA-3.0, via Wikimedia Commons)

If you imagine an extreme case, where your tRNA fragment is exactly divided in half (this is almost the case for your example!), or you by chance get only a large part of the variable loop, then you
see that you loose a large portion of the most significant base pairing (acceptor stem, anticodon loop+ stem). That will possibly drop the score immediately,
because the structure left would be a single hairpin loop, which are quite frequent and not significant for a tRNA. That means that for truncated sequences, blast might be more sensitive than RFam or tRNAscan-SE.
Question: I have a list of genes. and i want to extract the snps and indels from my VCF file (that i generated using GATK pipeline ) from genes coordinates on . The list of genes coordinates:

Gene Name Accession_no. Start_Position End_Position Strand
Rv0194 NC_000962.3 226878	230462 +

I was looking bedtools but it is asking for .bed format of genes nd as well .bed of bam files. how to do it ? or any other options/tools/scripts?

Like i tried tabix:

bgzip ERR038736_UnifiedGenotyper_variants_raw_snp.vcf
tabix ERR038736_UnifiedGenotyper_variants_raw_snp.vcf.gz
tabix ERR038736_UnifiedGenotyper_variants_raw_snp.vcf.gz AL123456.3:226878-230462 > Rv0194

and this gave me the variants like this:

AL123456.3 227098 . T C 6730.77 . AC=2;AF=1.00;AN=2;DP=172;Dels=0.
AL123456.3 228069 . G A 7132.77 . AC=2;AF=1.00;AN=2;BaseQRankSum=-
AL123456.3 228168 . G C 6682.77 . AC=2;AF=1.00;AN=2;DP=171;Dels=0.

But this is not a vcf file and i can only extract it one at a time. I want to extract all variants against a list of coordinates and store it in a vcf output.

Can anyone help me it this?

Answer: BEDOPS:

$ vcf2bed --insertions < snps.vcf > insertions.bed
$ vcf2bed --deletions < snps.vcf > deletions.bed
$ vcf2bed --snvs < snps.vcf > snvs.bed

https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html
Question: Hi, I have one BAM file which contains all alignments (include those not used in variant calling, such as non-PF, non-mapping and duplicate reads) generated for an assembly. How to filter these useless mapping? I know that Picard MarkDuplicates can be used to remove duplicates.Thank you.

Answer: You can use samtools to do this. e.g. to remove reads that did not align, you can do:

samtools view -F 0x04 -b in.bam > out.aligned.bam

to only include paired reads, use:

-f 0x02

Check the other bitwise flags on this page.

But you might not want to exclude those as they could be used for finding structural variations.
Question: I have a 20 nucleotide pattern. Now I want to find ALL of its occurences in a fasta file with about 200 sequences.

Then I want to ouput the positions either to a gff or into a sam file to be able to load into IGV. What is the most efficient way to do that?

It can have about 1-2 indels and about 2-4 SNPs. It is just important to find where the approximate pattern occurs often...

Answer: How about good old blat or gmap?
Question: Hi,

Is there any way to force the tree to look like a cladogram in ETE?

I can make the tree I want , but now I just want to add the style of the tree to have aligned branches the terminate at the same point.

Thanks

Answer: You need to convert the tree into ultrametric. ETE allows to do this with the `tree.convert_to_ultrametric()` method, which balances all branches in the tree to meet a given tree length. However, you can use a very simple strategy to add an extra offset to only leaf nodes so the rest of the tree conserves original branch distances.

from ete3 import Tree

#create a random tree
t = Tree()
t.populate(25, random_branches=True)

# prevent using labels and symbols in internal nodes
for n in t.traverse():
n.img_style["size"] = 0
n.img_style["vt_line_width"] = 0
t.render("regular.png")

#Convert to ultrametric
most_distant_leaf, tree_length = t.get_farthest_leaf()
current_dist = 0
for postorder, node in t.iter_prepostorder():
if postorder:
current_dist -= node.dist
else:
if node.is_leaf():
node.dist += tree_length - (current_dist + node.dist)
elif node.up: # node is internal
current_dist += node.dist
t.render("ultrametric.png")
Question: I am developing an R package and I plan to eventually upload it to BioConductor. However, after I went over the guidelines, I realized that it seems to be more restricted comparing with the requirement from CRAN. I followed the book: "[R Package][1]" by Hadley Wickham, but I believe this book is designed to help people submit CRAN packages, not too specifically for BioConductor.

Therefore I wonder if there is any special issue I should pay additional attention here. Especially formats I guess.

Thank you.

[1]: http://r-pkgs.had.co.nz/

Answer: In addition to all the good reasons and other issues mentioned in the other answers, I would point out to two resources that have not been explicitly mentioned so far. One is the [developers corner][1] in the bioconductor site, and inside it, the [package guidelines][2].

The developers site centralizes all the information regarding package development and includes the [release schedule][3], access to [build reports][4], the [new package primer][5] and much more.

In the guidelines many issues regarding the development of Bioconductor packages are discussed in detail. There is special mention about the issue of having [packages both in CRAN and in Bioconductor][6] (something not recommended). Also, it is advised that authors having questions about Bioconductor package development ask them in the Bioconductor development list ([bioc-devel][7]).

[1]: http://bioconductor.org/developers/
[2]: http://bioconductor.org/developers/package-guidelines/
[3]: http://bioconductor.org/developers/release-schedule/
[4]: http://bioconductor.org/checkResults/
[5]: http://bioconductor.org/developers/how-to/buildingPackagesForBioc/
[6]: http://bioconductor.org/developers/package-guidelines/#duplications
[7]: http://bioconductor.org/help/mailing-list/
Question: Hi

I recently asked a question as to whether exons can have different reading frames in the same gene. The answer was, surprisingly to me, yes. In this question I am assuming for simplicity that an exon has the same reading frame in all transcripts it features in.

Lets imagine a gene running from bases 1 to 200. Exon 1 is from bases 50 to 100 and exon 2 is from 125-150 and exon 3 is from 175-200.

Then lets say there is a splice variant (variant 1) containing exons 1 and 2 and another splice variant (variant 2) containing exons 1 and 3.

Then lets imagine there is a SNP in exon 1 at base 60. Could this SNP have a different effect in the different transcripts? For example is it possible that in variant 1 the exon is not translated at base 60 so the SNP is an an upstream untranslated region, but in variant 2 the SNP is in the translated region and causes a synonymous (or non-syn) amino acid change?

So another way of asking the question is this: if an exon appears in multiple transcripts, does it always play the same role in those transcripts? E.g is it always translated or not translated, and if it is translated, does translation start at the same place. I am assuming the answer to this will be 'no' because if an exon can have multiple reading frames in the same gene then all sorts of other scenarios are possible.

If what I describe does happen, roughly how common or uncommon is it?

Many thanks in advance

Answer: I've piped the program I wrote for your previous question and described [here](http://plindenbaum.blogspot.com/2010/10/where-are-alternative-reading-frames-in.html) to cross the data with dbSNP/UCSC:

java BioStar3034 | head 1000 | tr ":-" "," |\
awk -F ',' 'BEGIN { printf("select * from snp130 where chrom=\"_\" ");} {printf(" OR (chrom=\"%s\" and chromStart=%d)\n",$1,int($2)-1);} END {printf(";\n");}' |\
mysql -h genome-mysql.cse.ucsc.edu -A -u genome -D hg18

One of the result is rs1152725 http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg18&position=chr1%3A1209324-1209324

If the alternate base is 'C' instead of ref='T'

- transcript 1: codon `CTC (L) -> CCC (Proline)` (non-synonymous)
- transcript 2: codon `TCT(Serine) -> TCC (Serine)` (synonymous)

so my answer is YES: it happens :-)
Question: How does the kmer length affect the total assembly size and the N50 statistic with single end-reads and why?

Answer: Very generally spoken:

With larger kmer size there is a better chance of avoiding ambiguities in the graph between similar regions (repeats, paralogs,...). Ambiguities occur if kmer exist multiple times within the genome. (Unresolvable) ambiguities terminate contigs, hence larger kmer sizes in theory increase N50. However, large kmer sizes are much more sensitive to sequencing errors, heterozygosity and coverage.

Assembly size depends on how the assembler handles small ambiguities (bubbles) in the graph and how it handles low coverage paths. With small kmer sizes and some sensitivity to bubbles, you are more likely to generate a single contig for a slightly noisy region. This might be good in case of SNPs or bad if you merge repeats, .... This effect leads to smaller assembly size for small kmers. With large kmers you are more likely to generate different fragments for exiting variants, which increases assembly size. However, if these fragments fail for example internal coverage or length cutoffs, the final assembly may actually be smaller.
Question: Hey I have a fasta file which contains more than 80 sequences. I want to take each sequences in different files and the file name should be header of corresponding sequence.

Answer: http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/seqretsplit.html

seqretsplit
Question: Hi,

Being new to python scripting, I am facing a simple problem.

I am writing a python script "MyScript.py" in which I am calling the third-part program PyMOL through the import function.

Here is the simplified script (mainly inspired from this blog: http://doeidoei.wordpress.com/2009/02/11/pymol-api-simple-example/)

```
#! /usr/bin/python
# Usage: MyScript.py

# launch PyMOL from the terminal
import __main__
__main__.pymol_argv = [ 'pymol', '-qc']

import sys
import pymol

pymol.finish_launching()

# Input
file = sys.argv[1]

#run PyMOL command
pymol.cmd.do("load %s" % file)

#Exit PyMOL
pymol.cmd.quit()
```

I get this error message:

```
Traceback (most recent call last):
File "./MyScript.py", line 10, in
import pymol
ImportError: No module named pymol
```

My path to Python site-packages: `/Library/Python/2.7/site-packages`

My path to PyMOL: `/Users/me/my_apps/MacPyMOL_1.3.app/Contents/MacOS/MacPyMOL`

(I made an alias `pymol='MacPyMOL'` in `.bash_profile`)

I also tried to replace in the script:

```
import sys
import pymol
```

by

```
import sys
sys.path.append("/Users/me/my_apps/MacPyMOL_1.3.app/Contents/MacOS")
import MacPyMOL
```

but still the same.

I looked at other forums but it didn't solve the problem:

- https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=CCP4BB;24c8090b.1011
- https://www.mail-archive.com/pymol-users@lists.sourceforge.net/msg08176.html

If someone could direct me to a solution...

Thanks!

Answer: Hi,

In MacOSX with a default installation of MacPymol you cannot use the API because there is no _cmd.so file. What I usually do it to install pymol via fink/macports or compile it from source somewhere. Not the best solution but then you can get the file you are missing.

If you manage to find a better solution, it'd be great to share !