cannin/biostars_qa huggingface dataset:
Question: Hello all,

I've read several papers estimating the quality of ChIP-seq and most estimate percentage of peaks that include the consensus sequence at about 50%. Having lower FDR does not improve this very dramatically (to 70% maybe).

So what are the high confidence peaks that do not contain the target motif? What causes them? Is there a systematic study of any sort?

I'd be grateful for any pointers.

cheers

Answer: Some possible explanations could be:

Secondary binding/protein complexes: Another protein is binding the DNA and your ChIP:ed protein is binding to that protein.

Generally open chromatin: an active regulatory region (perhaps "unlocked" by a "pioneering TF") is easier to bind to in the absence of a motif.

Degenerate motif: there is actually a binding motif but it is relatively dissimilar from the consensus

Artefactual regions, see e g the "black list" here: https://sites.google.com/site/anshulkundaje/projects/blacklists

Since you have an FDR it seems you have a control library so the following may not apply to your case but is interesting in general: even in negative control libraries, regions like transcription start sites give robust peaks: http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0005241
Question: Hi there!

Can anyone help in answering how to split single .txt file with multiple fasta sequences into multiple files containing each of the fasta sequence?

Thank you for your time and help!

Answer: using awk

awk '/^>/ {if(x>0) close(outname); x++; outname=sprintf("_%d.fa",x); print > outname;next;} {if(x>0) print >> outname;}' *.fasta
Question: Hi! I'm new to bioinformatics, but am working with some fastq files that have some strange base quality distributions, see image. It is strange as we see only 4 unique phred scores across the whole file, which seems surprising given illumina sequencing has 41 possible scores. This is happening across multiple files, and these files are straight from the sequencing company. The values correspond to phred scores of "F" "," ":" and "#".

I have confirmed this behaviour with multiple people in my team, so this is not an analysis problem this is an issue with the files (also obvious when looking at raw read phred scores).

I also found [this][1] other question on Biostars, which it's hard to tell but they seem to have the same behaviour, suggesting perhaps it is a common issue. Does anyone have any idea what is happening? The reads themselves seem normal when compared to reference genome.

We have contacted the sequencing company and they haven't really provided clarity so thought that maybe people here could provide some insight.

Thanks in advance!

![Aggregation of base qualities over fastq files for same sample][2]

[1]: https://forum.qiime2.org/t/fastq-quality-too-good-to-be-true/12674
[2]: /media/images/8e68fa10-5a2a-48e8-a38e-1872eb33

Answer: Reposting from @lievensterck comment:

"If I remember correctly, was Illumina not gonna change it's qual scores, in binned approach to reduce file size?"

From which I found [this][1] answer which seems to be it.

[1]: https://www.biostars.org/p/408032/
Question: I want to uninstall *Bioconductor* [affycoretools package](https://www.bioconductor.org/packages/release/bioc/html/affycoretools.html) and reinstall it, who can help me?

Answer: Remove a package with `remove.packages()` e.g.

remove.packages("affycoretools")

**Affycoretools** is a Bioconductor pacakge, so reinstallation needs their install script / the **BiocInstaller** package e.g.

source("https://bioconductor.org/biocLite.R")
biocLite("affycoretools")

---------

**Edit:** To install this package, start R (version "4.0") and enter:

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")

BiocManager::install("affycoretools")
Question: Hi, I'm an undergraduate student. Please help me.

I want mapping-bam-file from sra-dataset from NCBI in order to analyze heterogeneity of mouse ESC.

The dataset is generated from a paper below.(GSE60749)

*Roshan M.Kumar et al. Deconstructing transcriptional heterogeneity in pluripotent stem cells. Nature(2014)*

Please teach me how to get adapter sequence used in this experimentation

and what I should use for quality control(Prinseq?ShortRead?).

For example, I want to try GEO Sample [GSM1486817][1] sra file.

Can anyone give me process of quality control of this sra?

I thank you for reading it through.

Any help will be appreciated.

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1486817

Answer: You can find the adapters used by reading a bit and googling around. "Nextera XT DNA Sample preparation reagents (Illumina)" were used to prepare the samples (as found here under "study summary", or here under "library").

You can check for adapter contamination with FastQC.

edit: if the downstream analysis to be performed is mapping with BWA, Bowtie2 or any mapper which performs local alignment, adapter should not have a major impact. Besides, the reads he pointed to in his question are 25+25, which should be shorter than Nextera typical insert size, and adapter contamination should be really low.
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: http://manuals.bioinformatics.ucr.edu/home/R_BioCondManual#TOC-Vectors

http://www.bioinformatics.babraham.ac.uk/training.html#rintro
Question: Dear all,

Does STAR 2.6 `--quantMode GeneCounts` only quant unique mapping reads? How to specialize it only quant the unique mapped reads?

Thanks a lot.

Answer: Yes, and you will have to choose a column among the results, please refer to https://www.biostars.org/p/218995/
Question: I have a dataframe in R that looks like this:

V1 T1 T2 T3 T4 T5
CXCL6 0.8536601 1.0903336 3.7633042 5.5800459 5.8477150
PPBP 0.7739450 0.3587961 0.5073359 0.2743522 0.6221722
CXCL10 0.1258370 -0.3535165 -0.7460387 3.5604672 0.1971432
CXCL11 -0.2563139 0.7117200 0.0000000 -0.2288303 0.9955557
CXCL12 0.6181279 1.7529310 1.7637760 1.2752787 1.2284810

I want to keep the rows that have values only between -1 and 1.

I have tried this command but unfortunately it does not work.

condition1 <- Genes[,c(2:6)] <=-1 & Genes[,c(2:6)] >=1
Genes <- Genes[condition1,]

Can someone tell me where I am wrong so that I can successfully filter my dataframe.

Answer: There should be more tidy way for that, but here is my solution.

rowIdx = as.numeric()
for(i in 1:nrow(Genes)){
cond = Genes[-1][i,] >= -1 & Genes[-1][i,] <=1
if(all(cond == TRUE)){
rowIdx = c(rowIdx, i)
}
}

Then subset your dataframe by `rowIdx`:

GenesSubset = Genes[rowIdx,]
Question: Dear all,

i have this table. does someone know how to de-aggregate it ?

plant	tissue	count
tomato	leaf	1
tomato	root	4
tomato	shoot	5
solanus	leaf	3
solanus	root	2
solanus	shoot	4

what i want is a dataframe like that:

leaf	root	shoot
tomato	1	4	5
solanus	3	2	4

thanks in advance for you tips

Answer: Buenas tardes amiga/o,

You can use `dcast()`:

df
plant tissue count
1 tomato leaf 1
2 tomato root 4
3 tomato shoot 5
4 solanus leaf 3
5 solanus root 2
6 solanus shoot 4

require(data.table)

dcast(data = df, formula = plant ~ tissue, value.var = 'count')
plant leaf root shoot
1 solanus 3 2 4
2 tomato 1 4 5

Kevin
Question: Hi~
I'm currently using 10x cellranger to analyse single cell RNA-seq data. According to their algorithm, reads mapping confidently to more than one exons will be discarded. However, there are paralogous genes in the genome that are largely identical and all the reads for such genes are discarded. Therefore, I was wondering if there is a way to change the algorithm to count the first (or a random) confident alignment. Unfortunately, I wasn't able to locate the file containing the algorithm. Any hints would be appreciated. Thanks every one!

Answer: The cellranger UMI deduplication algorithm does not handle reads that map among multiple genes, there is no "easy" way to handle this situation. You may be interested in taking a look at our quantification tool, [alevin](https://www.biorxiv.org/content/early/2018/10/24/335000), which we've designed, in part, to help deal with these cases. In addition to having a methodology for handling reads that map between multiple genes, it is much faster than cellranger.
Question: 0
down vote
favorite
I was experimenting Prokka and RAST annotation tools. So, I took a well-annotated swinepox virus genome from genebank *(NCBI Reference Sequence: NC_003389.1)*.

I ran those sequences on Prokka and RAST Seed server at the same time. I can see that only a few **(may be around 1%)** of the genes were annotated. Most of them were predicted as **hypothetical protein**. And the results were comparable between Prokka and RAST.

I would assume that these tools look for similar sequences in NCBI and find the best-match protein. But looks like that is not the case. It should be able to find that well annotated swinepox virus genome in genebank and predict most of the proteins.

Also, if almost all the genes are predicted as hypothetical protein then there is not much difference between gene prediction tool like Genemark and genome annotation tool. Are there any better annotation tools ? Or this is what we get ? Or have I misunderstood the concept of annotation ? Please someone help me to understand this thing.

*I have attached the image for comparison. Left one the Swinepox genome in .gb format and right one is the same genome annotated with Prokka.*

Image link : https://github.com/lrjoshi/sample/blob/master/docs/annotation.PNG

![Image][1]

[1]: https://github.com/lrjoshi/sample/blob/master/docs/annotation.PNG

Answer: > Also, if almost all the genes are predicted as hypothetical protein then there is not much difference between gene prediction tool like Genemark and genome annotation tool.

Genome annotation is a two-step process, first you have to predict where the genes are - which tools like Augustus and GeneMark do - then you have to assign function to the predicted genes - usually by means of similarity searches using good quality databases. This is what Prokka and RAST are doing, but integrated in a pipeline.

I don't know how RAST works, but prokka uses several programs to predict genes (protein coding, non-coding RNA, tRNA, rRNA, and more) from the genome, then, after having these predictions, it tries to annotated them by searching available databases. Prokka uses Prodigal for protein coding gene prediction, which I don't know if it appropriate for virus gene prediction. But the most important reason your annotation came up mostly as hypothetical proteins probably is you don't have installed an appropriate database to annotate viral genomes - you can pass one at run time with the `--proteins` option, which will have precedence over other installed databases.
Question: I have genotype data in the form of a number of birdseed files (one for each sample), of the form:

```
Composite Element REF Call Confidence
SNP_A-2131660 2 0.0060
SNP_A-1967418 2 0.0281
SNP_A-1969580 2 0.0074
SNP_A-4263484 2 0.0104
SNP_A-1978185 0 0.0034
```

I've loaded these into a matrix in R to perform QC upon them and annotate them with the rs ids instead of the probe id, resulting in something like this:

```
> geno[1:4,1:4]
TCGA.A7.A0D9 TCGA.A7.A0DB TCGA.A7.A13G TCGA.AC.A2FB
rs2887286 2 2 2 2
rs1496555 2 2 2 1
rs3890745 2 2 1 2
rs10489588 1 0 0 0
>
> dim(geno)
[1] 693085 94
```

I want to use 1000 Genomes data to impute genome wide SNP data, however I'm not sure how to get from this format to the format required by IMPUTE2. I've been reading through [this wiki][1] and it seems like it might be easiest to convert this to PED format and then use GTOOL to convert to IMPUTE format.

Are there any existing tools that can do this reformatting directly? Or will I have to write a script to do the matrix > PED reformatting and then use GTOOL for PED > IMPUTE?

[1]: http://genome.sph.umich.edu/wiki/IMPUTE2:_1000_Genomes_Imputation_Cookbook

Answer: If you use shapeit2 to prephase the chromosomes (best practice) it will accept bed files. Shapeit2 outputs files that can be read by Impute 2. And shapeit2 can read bed files directly, with the `--input-bed` flag.

The filetype impute2 uses is called Oxford, and with plink1.9 you can convert most filetypes to Oxford format. I've never heard of birdseed format, but here is a script that can convert birdseed to bed: https://www.broadinstitute.org/mpg/birdsuite/howto.html
Question: Hi guys,

Question, I have a couple of sam files that I got from aligning with Bowtie2. I am trying to sort through them to remove the usual stuff (multi-mapping etc etc) after reading posts and guides etc etc. So far I've been successful in taking my sam files, running `sed '/XS:/d' algn.sam > mappedonce.sam` and getting rid of all those with an XS flag. When I then run the following :-

```
samtools view -c -f 4 mappedonce.sam
749598

samtools view -c -F 4 mappedonce.sam
12116420
```

I can see that I have 749598 reads that are unaligned which I want to remove. So I ran this as per posts on biostars

samtools view -F 4 -o mapped1_mapped.sam mappedonce.sam

However, when I then take this file after running the -F 4 filter, if I try and run anything else such as flagstat or even the `samtools view -c -F 4 mappedonce.sam` again, I get the following error:

```
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 1
[main_samview] truncated file.
```

So what is samtools doing to my file? Is there a reason why samtools would remove a header line and truncate my file? Am I doing something wrong? Consequently, if anyone has any non-samtools solutions to remove these 749598 reads I'd be open to that!

Answer: You need the `-h` option in `samtools view` so the header is also printed out.
Question: Hi, all.

I will output MDS plot with ggplot2.

But, I couldn’t read the data because the class of the data wasn’t support in ggplot2.
Is there a way to convert “dist” to “data.frame”?

Although I tried “as.data.frame()“, I was not able to do.
The following error messages was outputed.

> as.data.frame(d)
Error in as.data.frame.default(d) :
cannot coerce class ‘"dist"’ to a data.frame

I'm a beginner in bioinformatics, so sorry for the basic question.
Could you please advise me?
Thank you very much for your appreciated!

Answer: You have a specific function designed to perform this in the `metagMisc` package :

devtools::install_github("vmikk/metagMisc")
library("metagMisc")
dist2list(d, tri = TRUE)

You have also plenty other solutions as you may find [here][1].

[1]: https://stackoverflow.com/questions/23474729/convert-object-of-class-dist-into-data-frame-in-r
Question: I am thankful for Obi Griffith sharing his code on how to make heatmaps in R where you can add row and column information together with the color mapping information. I use it all the time to make pictures of hierarchically clustered expression profiles measured by expression arrays.

Usually I scale the input matrix on rows (probesets) to Z-scores in order to make weigh equally in the clustering. The colors of the heatmap are mapped linearly to the Z-scores. \

The problem is, when you have a few extremes in your data set, they tend to make your heatmap a little faint. In the attached heatmap you can see that the Z-scores range from somewhere around -10 to 10, however, about 99% is in the range of -4 to 4. I would like to change the way colors are mapped to the Z-scores in such a way that the green-to-red goes from -4 to 4 and the colors of the Z-scores outside this range just do not become any more green or red then they are already, any ideas?

# import gplots for greenred colors
library(gplots)
# generate random matrix
x<-replicate(100,rnorm(100,0,1))
# draw heatmap
heatmap.3(x,density.info="histogram",col=greenred)
# introduce outliers
x[1,1] <- -10
x[2,2] <- 10
# draw heatmap with outliers
heatmap.3(x,density.info="histogram",col=greenred)

Answer: You can use the breaks argument to do this, from the heatmap.2 documentation:

breaks
(optional) Either a numeric vector indicating the splitting points for binning x into colors, or a integer number of break points to be used, in which case the break points will be spaced equally between min(x) and max(x)

To use to achieve your desired effect:

#breaks for the core of the distribution
breaks=seq(-4, 4, by=0.2) #41 values
#now add outliers
breaks=append(breaks, 10)
breaks=append(breaks, -10, 0)
#create colour panel with length(breaks)-1 colours
mycol <- colorpanel(n=length(breaks)-1,low="green",mid="black",high="red")
#now plot heatmap - I'm using heatmap.2, but am assuming this will work with Obi's heatmap.3
heatmap.2(x,density.info="histogram",col=mycol, trace='none', breaks=breaks)

The major drawback is that the outliers are a lot less obvious, but clearly the heatmap overall is a lot less 'washed out'.

Results:

Before

After
Question: Hi BioStars,

I am using pysam to count read from bam file. I use pysam.AlignmentFile to read bam file and pysam.AlignmentFile.count to do read count, but it requires bai file which is an index should in the same path where bam file is. However I have no writing permission in bam file path. Is it possible to input index bai file from other path, I mean save bam file and its bai file in different path? Thanks very much if some guys can help this!

Answer: Easiest is probably to create a symbolic link to the bam (`ln - s`) in a directory where you can write, or `/tmp/`, and write the index there. You can probably fix this using the `os` module
Question: Hello all,

I have barcode count data corresponding to the viability of 25+ pooled bacterial strains under various conditions. The marginal distribution of untreated strain counts appears to be Negative Binomial.

I'm trying to use DESeq2 to analyze these data, using a matrix of strains ("genes") as rows and conditions as columns. Since the variation of counts between most conditions for most strains is very large, but between replicates is relatively small, it seems sensible to estimate dispersions (in this case) on a gene- and condition-wise basis.

The language in the DESeq2 vignettes and pre-print seems to suggest the dispersion estimates are "gene-wise". So if you run `DESeq()` followed by `plotDispEsts()`, each point corresponds to the variance estimate of a gene across conditions (in my case, strain), or the variance estimate between replicates of a gene under one condition?

I think the conceptual difference I'm talking about is the same as that between `blind=TRUE` and `blind=FALSE` in the `rlog()` and `varianceStabilizingTransformation()` functions.

Finally, if DESeq2 does estimate dispersions on a solely gene-wise basis, would it be reasonable for me to estimate the dispersions of my data subsetted by each condition in turn, and then feed those results into my whole `DESeqDataSet` object using `dispersions()`?

Many thanks for taking the time to read, and for any suggestions you might have.

Eachan

Answer: The *estimateDispersion* function will estimate gene wise dispersions using all columns. Although a GLM fit (based on within group variances and means) is done before estimating so the differences in conditions are also incorporated. So I think you don't have to (and probably shouldn't) use condition-wise estimation.

By the way, why do you want to use DESeq2 here? I mean with a dataset with 25 rows only..
Question: Are there any publicly available raw sequencing data for COVID-19? In BAM for FASTQ format?

Answer: `sra-explorer` [search with][1] `COVID`. Download using instructions here : https://www.biostars.org/p/366721/

[NCBI's SARS-CoV-2 genome's page][2] has direct links to SRA accessions (scroll down the page).

[1]: https://sra-explorer.info/
[2]: https://www.ncbi.nlm.nih.gov/genbank/sars-cov-2-seqs/
Question: Hi Everyone,

I recently heard about the genome in a bottle project.

I found raw fastq data for NA12878 samples from their GitHub repository. Now I'm trying to reproduce their analysis.

https://github.com/genome-in-a-bottle/giab_data_indexes/blob/master/NA12878/sequence.index.NA12878_Illumina_HiSeq_Exome_Garvan_fastq_09252015

I would like to know if it's available a reference code to compare my pipeline (not only the final result but the code).

Thank you in advance,

AY

Answer: Hi Yussab,
as you did not replay, so i do not know what are you planning to do exactly.
it seems that you do have your own pipeline which mean you have the final VCF where it's generated from NA12878.
i know these 2 tools to compare you result with GIAB Data set. (I do not know if there is something else available)

1- Hap.py (https://github.com/Illumina/hap.py):

hap.py truth.vcf query.vcf -f confident.bed -o output_prefix -r reference.fa

2-RTG tools (https://github.com/RealTimeGenomics/rtg-tools):

generate SDF file from Fasta file (Your reference):

rtg format "Your_Fasta_File".fasta -o "File_Name_ToOut_Put_Result"

Compare the tow VCF file:

rtg vcfeval -b truth.vcf -c query.vcf -o Directary_To_OUTPUT_THE_RESULT -t PATH_TO_SDF_FILE -e GIAB_BED_FILE.bed --region= Your_Bed_File.bed

I highly recommend you to read the tools manual for better understanding.
Question: Hello everybody

Can any one help me understand the MDS plot for two groups of samples like control vs. treated we generate from edgeR or Deseq2. I understand the figure as it separates the control from treated very well, but if we have the y-axis as the logFC, I don't get what the x-axis is? What is the x-axis in such 2D figure?

Thanks

Answer: Hello, the y-axis in a MDS plot does *not* represent logFC, and neither does the x-axis.

Provided you have generated the plot according to the workflow, RNA-seq workflow: gene-level exploratory analysis and differential expression, on Bioconductor, your x- and y-axes will be representative of **Euclidean distances** between your samples. These Euclidean distances will be produced/repeated across multiple dimensions in different ways, but the standard way of representing MDS is to just plot the Euclidean distances from the first 2 dimensions, with x-axis being *Dimension 1* and y-axis being *Dimension 2*. The dimensions are ordered based on how well they fit your samples (for further reading, take a look at the cmdscale *Details* section.

MDS is commonly used in genetic studies to find relationships between samples based on genotype, but, as you can see, it's also used in RNA-seq. For information on a related method, i.e., principal components analysis, please see my threads here:

- https://www.biostars.org/p/280615/#280634
- https://www.biostars.org/p/282685/#282691
- https://www.biostars.org/p/271694/

-------------------

For what it's worth:

- A volcano plot represents logFC (x-axis) and negative log10 P- or
adjusted P-value (y-axis)
- A MA plot represents average expression (x-axis) and logFc (y-axis)