cannin/biostars_qa huggingface dataset:
Question: Dear users,
I read through UCSC wiggle format (BigWig) in depth. I understand the format about variable step, fixed step and use of bigwig to visualize in a browser etc. However, I don't have clarity on using ChIP-Seq, ATAC-Seq data in bigWig format.

> Example of the Wig file after converting from bigwig format using bigWigToWig. This example data is from ATAC-Seq.

#bedGraph section chr1:0-870999
chr1 0 9999 0

chr1 9999 10099 16.561

chr1 10099 10199 24.2045

chr1 10199 10299 2.54784

chr1 10299 10399 5.09568

chr1 10399 10499 11.4653

chr1 10499 10599 7.64352

chr1 10599 10699 3.82176

chr1 10699 13199 0

In this context, could someone help understanding:
1. What is the real number value in 4th column. Is the read depth for that position? or some transformed value of read depth. Typically in ChIP-Seq or ATAC-Seq, what is the value that one would represent in this column.

2. If this is threshold, how users specify thresholds for selecting significantly enriched regions. (I am not sure if any statistical test is associated with significance, I cannot find any reference but users call it so).

3. I have 6 ATAC-Seq bigwig files for 6 different samples. How do I find the regions of interest in at least 4 of samples.

Thank you for your help.

-Adrian

Answer: The bedGraph (or bigwig) format is always the same: `chr-start-end-value`. Value can actually be anything that can be associated with a stretch of DNA as defined in the first three columns. It can be the raw read count for that interval, it can be normalized read count like reads per million, it can be an enrichment score for this experimental condition over a control experiment, it can be the mean methylation store, the GC content etcetc. Most commonly, people use bedGraph/bigwig to create browser tracks displaying the normalized read count across the genome, and in this case, it would not matter if it is ATAC-seq, ChIP-seq or Whatever-seq. One simply counts the number of reads that cover each base and aggregates bases with equal coverage into one interval to make the files smaller, so if the first 100 bases of a chromosome have coverage of 0, one would write:

chr1	0	100	0

instead of 100 intervals like:

chr1	0	1	0

chr1	1	2	0

chr1	2	3	0

(...)

For statistical analysis, one typically calls peaks (e.g. with MACS) and then makes a count matrix to obtain the raw counts for each replicate per peak. Significances between conditions are then inferred with appropriate statistical frameworks, such as `DESeq2, edgeR, csaw` etc. Please use the search function and google for differential analysis of ATAC-seq data, there is plenty of material available.
Question: Hi,

The fastq format right now has: header, sequence string, "phantom" header and quality string. For storage purposes, why doesn't the fastq format incorporate nucleotide information in the quality strings? Is it just to make it more human readable?

Answer: FASTQ is almost never stored uncompressed. With compression, merging bases and quality may actually hurt compression as compression is usually less efficient when you mix different types of information together. FASTQ is not meant to be read by a human from the start to end, but it is meant to be eye-read in a small portion and manipulated with the many unix tools. These are critical features.
Question: "discontiguous megablast" on the NCBI website, as well as other implementations, appears to have a "Discontiguous Word Options" parameter set. These include a template length and template type; the latter can be:

Coding

Maximal

Two templates

What do these options do?

Answer: The way BLAST works is by first matching a word from the query and the database, this match is then expanded using dynamic programming. The word is usually 7-256 in length. When using discontiguous BLAST the word is 11-12 letters in a 25 letters long sequence. The way these letters that must match between the query and database are spread differently along the longer sequence. In coding the matching letters will be in the first two nucleotides of every triplet, in maximal they will be spread in a predefined pattern that should maximize the number of matches, two templates just try to match these two patterns.

In short discontiguous BLAST allows matches that don't have a 7 letters word that perfectly match between the two.
Question: I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

```
>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence
```

And I have a list like the following:

```
Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11
```

I want to retrieve those sequences that match with names present in the list.

Answer: You can do this easily by installing the FAST: Fast Analysis of Sequences Toolbox ([publication][1])([github][2]).

You can install it by using this command (only use sudo if you need to):

(sudo) perl -MCPAN -e 'install FAST'

Once it is installed here is a small bash script to do what you need:

```
#!/bin/bash

while read line ;do
cat original_fasta.fa | fasgrep -di "$line" >> reduced_fasta_file.fa
done < species.txt
```

[1]: http://journal.frontiersin.org/article/10.3389/fgene.2015.00172/abstract
[2]: https://github.com/tlawrence3/FAST
Question: Hi friends,

I am attempting to sort my bam files that I obtained from my bowtie sam files. I am not indexing them appropriate according to this error I am receiving after creating my bam file.

random alignment retrieval only works for indexed BAM or CRAM files.

I understand I am suppose to index the file before sorting them.

#creating the appropriate files
samtools view -Sb sample.sam.pair > sample.pair
samtools view -bt ~/bigdata/refgenome/genome.fa.fai - - | samtools sort sample.pair -o sample.pair.bam

samtools view -Sb sample.sam.single > sample.single
samtools view -bt ~/bigdata/refgenome/genome.fa.fai - - | samtools sort sample.single -o sample.single.bam

#merge
samtools merge sample.all.bam sample.pair.bam sample.single.bam -@ 2
rm sample.pair sample.single

#index the final bam
samtools index sample.all.bam

Any help would be appreciated.

Answer: I think you're over-thinking things :)

You can only index BAM files on position, and only when the data is sorted by position to begin with (don't ask...)
So to sort by position just do:

samtools sort my.sam > my_sorted.bam

Then index with

samtools index my_sorted.bam

It's as easy as that. If you want to merge the output files from bowtie do that as the very first step, because I don't think samtools performs any optimisations for merging sorted BAMs/SAMs. However, i'd also recommend against bowtie2 in favour of STAR or BWA-MEM, but that's just a personal preference at the end of the day.
Question: Hello,

I've hit a bit of a snag here: I had to update to the latest version of R to get ballgown working, but one of the functions I see included in filtering out low expression genes, rowVars, is not available for this version of R.

Using this tutorial: https://rpubs.com/kapeelc12/Ballgown

The rowVars function is employed in filtering:

bg_filt = subset(bg,"rowVars(texpr(bg)) >1",genomesubset=TRUE)

I have R studio, would it be possible to download an older version of R and switch to it temporarily to use rowVars, then switch back to continue the ballgown output analysis?

If so does anyone know what the last version of R rowVars worked on?

Or should I try to find an older version of ballgown, and downgrade my R version? I know there are many functions on bioconductor that do not necessarily work on the latest version of R, so is there a benifit to running older versions of R generally when working with these packages?

Answer: Is it the `rowVars` function from [metaMA][1] ?

install.packages('metaMA')
library(metaMA)
?rowVars

[1]: https://www.rdocumentation.org/packages/metaMA/versions/3.1.2
Question: I have 2 sorted Bam files and I want to merge them and get a VCF file from them using samtools pileup.

first i merge these 2 files with `samtools merge` and then I get a VCF file as output.

samtools mpileup genome.fa merged-bam.bam output.vcf

then I use

samtools mpileup genome.fa 1.bam 2.bam output.vcf

but the output of these 2 approach had different size. where am I doing wrong?

thanks a lot

Answer: my answer depends of your sam header+ read-group ( https://gatkforums.broadinstitute.org/gatk/discussion/6472/read-groups )

1) samtools merge : would produce one bam and you'll get only **one** virtual sample in the final vcf (see the columns after FORMAT in the '#CHROM' line)

2) the 2nd command 'samtools mpileup genome.fa 1.bam 2.bam output.vcf' will produce a vcf with **two** samples (this is what your want in 99% of the use cases)
Question: Dear all,

considering a RNA-seq experiment and analysis that provides the expression values as TPM, please would you let me know what is a minimum TPM value in order to consider a gene to be expressed ?

talking about RPKM.FPKM units, I remember that a gene was considered expressed if RPKM (or FPKM) > 1 ... thanks a lot,

-- bogdan

Answer: As already pointed out, there is no ideal cutoff. However, there is at least one method, zFPKM, that tries to define an expression cutoff.

BioC: https://bioconductor.org/packages/release/bioc/vignettes/zFPKM/inst/doc/zFPKM.html

Publication: https://www.ncbi.nlm.nih.gov/pubmed/24215113

> the community adopted several heuristics for RNA-seq
> analysis, most notably an arbitrary expression threshold of 0.3 - 1
> FPKM for downstream analysis. However, advances in RNA-seq library
> preparation, sequencing technology, and informatic analysis have
> addressed many of the systemic sources of uncertainty and undermined
> the assumptions that drove the adoption of these heuristics. ... We
> use ENCODE data on chromatin state to show that ultralow-expression
> genes are predominantly associated with repressed chromatin; we
> provide a novel normalization metric, zFPKM, that identifies the
> threshold between active and background gene expression; and we show
> that this threshold is robust to experimental and analytical
> variations.
Question: Hi all, I'm trying to run a single-cell velocity analysis on hg19 data using kallisto, following along with this tutorial: [https://www.kallistobus.tools/velocity_index_tutorial.html][1]

In the tutorial, the introns file they generate theoretically looks like this, with the familiar looking ENST transcript IDs:

$ head -4 introns.bed
chr1	12118	12721	ENST00000456328.2_intron_0_109_chr1_12228_f	0	+
chr1	12612	13329	ENST00000456328.2_intron_1_109_chr1_12722_f	0	+
chr1	11948	12287	ENST00000450305.2_intron_0_109_chr1_12058_f	0	+
chr1	12118	12721	ENST00000450305.2_intron_1_109_chr1_12228_f	0	+

However, mine look like this, though I followed their exact instructions on how to generate the introns.bed.gz file from the UCSC table browser.

>> head -4 introns.bed
chr1	12200	12639	uc001aaa.3_intron_0_27_chr1_12228_f	0	+
chr1	12694	13247	uc001aaa.3_intron_1_27_chr1_12722_f	0	+
chr1	12200	12672	uc010nxr.1_intron_0_27_chr1_12228_f	0	+
chr1	12670	13247	uc010nxr.1_intron_1_27_chr1_12698_f	0	+

Because of this, I don't think I am able to continue following along with the tutorial, and also don't imagine my file is correct anyways. Can someone please help me understand a) why I can't get normal-looking transcript IDs for GRCh37 (hg19) from UCSC or b) whether I am just doing something totally wrong or c) whether this is fine and there is a workaround. I am just confused as to what those strange IDs are...

Thanks for the help!

[1]: https://www.kallistobus.tools/velocity_index_tutorial.html

Answer: If you want your transcript IDs to have the ENST identifiers, you should select Ensembl Genes instead of UCSC Genes for the track.
Question: What happened with miRecords mirna database website? The webpage is not available and I used it's data as part of an analysis in my paper. Now the reviewer can't see the page and asks for an explanation.

Edit:

Sorry, I forgot to post the link. The miRecords was suppose to be available at http://miRecords.umn.edu/miRecords.

Answer: The persistence (or lack of) bioinformatics resources online has been analysed a number of times:

- [On the persistence of supplementary resources in biomedical publications][1]
- [A survey of the availability of primary bioinformatics web resources][2]
- [404 not found: the stability and persistence of URLs published in MEDLINE][3]
- [Persistence and Availability of Web Services in Computational Biology][4]
- [A Survey of the Availability of Primary Bioinformatics Web Resources][5]

And sadly you have fallen foul of a resource that has either moved (with no redirect), been deprecated, or the server in the server room has died, and the one person who knew how to maintain the resource left the lab 3 years ago.

Reliance on sources that you cannot download and archive yourself, or are hosted outside of institutional repositories, or specific data repositories such as Figshare is asking for trouble. You will probably need to alter the paper to reflect data from a similar resource or attempt to contact the authors of the original paper. Or share some of the URLs above with your reviewer.

[1]: http://www.biomedcentral.com/1471-2105/7/260
[2]: http://www.sciencedirect.com/science/article/pii/S1672022907600175
[3]: http://bioinformatics.oxfordjournals.org/content/20/5/668.abstract
[4]: http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0024914
[5]: http://www.sciencedirect.com/science/article/pii/S1672022907600175
Question: I have a dataset of proteins that I have blasted against the uniprot-swissprot database.

I'd now like to identify which proteins are likely to have a mitochondrial sub-cellular localisation based on the sub-cellular localisation of their best blast hit in the swiss-prot database.

The fasta headers of the uniprot proteins look like this:

">sp|Q64602|AADAT_RAT Kynurenine/alpha-aminoadipate aminotransferase, mitochondrial OS=Rattus norvegicus GN=Aadat PE=1 SV=1"

I have found a gene ontology mapping file (link below) but the fasta headers don't contain the GO IDs necessary to map them.
ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/external2go/uniprotkb_sl2go

Is there some intermediate file that I need to use and does anyone know where to find it? Any help would be appreciated.

Answer: using xslt:

$ awk -F '|' '/^>/ {printf("%s\n",$2);}' input.fa | while read ACN ; do curl -s "https://www.uniprot.org/uniprot/${ACN}.xml"| xsltproc transform.xsl - ; done

Q64602 Mitochondrion

with transform.xsl:

https://gist.github.com/lindenb/92ae5d03183d1ff56a17684d30dd8f7e
Question: I have used tophat2 to map rna-seq reads to a draft genome. The alignment percentage is around 75-80% for all samples. When I take the unmapped reads and blast them, they hit the same organism, indicating the unmapped reads might have potential information. How do I deal with the unmapped reads and include them in DE analysis or any other downstream analysis ? Should I go with entirely different pipeline like trinity ?

Answer: I have tried STAR and the mapping percentage increased up to 90-92% ( with tophat2, it was only up to 75-85%). I will try BBMap soon.
Question: It seems that the multiple test for the GO enrichment result is important. However, the TopGO didn't do that, so I decide to do that by myself. The delemma looks like this: fdr correction with `p.adjust()` in R need the complete GO enrichment result, whereas the `GenTable(,topNodes =N)` in TopGO only give out the top N result. So how can I do the "fdr" mutiple test correction to the GO enrichment results from "TopGO"?

Sincerely

Answer: Instead of specifying `topNodes` to some value N, you can get all available GO Terms using:

```r
allGO = usedGO(object = GOdata)
# use it in GenTable as follows:
GenTable(GOdata, ... ,topNodes = length(allGO))
```

Then you can go about doing your p-value adjustment.
Question: Hi!
I need advice on the processing steps for my project. It would be really nice to get some feedback on this.
I want to perform Trinity de novo assembly from metatranscriptome samples. Since there are different microorganisms inside, I want to use the big assembly as a reference to map then the individual samples (replicates under 2 different conditions) and count the transcripts for testing differential expression.
The problem that I am struggling with is: for the big assembly I will have around or more than 200 million reads (PE), depending on how I will process the sequences (I could have 2 big assemblys, each one for the different conditions, and in this case will be less data; or I could maybe get a better "reference assembly" by using all data together) . So I don't now if it will be possible to perform this without problem on the requested resources using Trinity in a HPC cluster, until now I have only used 40 million as a maximum for assembly and its really difficult to keep running a job for so long.
Maybe you could give me some advice on how could I improve the data (pre)processing steps?

Another thing that I'm not sure about is: as I dont have a reference genome and I'm not expecting to have a big percentage of further annotation, It would be better for me to use merged PE (longer reads) that its about 20-30% of my sequences, but in this case I would loose the rest of the information from the unpaired reads AND I should treat my sequences in single mode with Trinity...
Is there a way that I could combine my merged data with the unmerged and include everything in my analysis? without having to treat everything as single end data?
Thanks in advance!

Answer: You should make only one consensus transcriptome with all your reads. After that, make the differential gene expression, both analyses can be made using Trinity. Overmore, you should use Trinotate for functional annotation. You will need to use a cluster with approx 250G of RAM during 24 hours. Follow the instructions in Trinity GitHub wiki: https://github.com/trinityrnaseq/trinityrnaseq/wiki.

Additionally, you should use Megan to analyze the putative species.

I have the pipeline to do that already so if you need any help contact me.
Question: Hello everyone,

I am interested in a sequence ranging from 111614 to 111868 in a fasta sequence (scaffold sequence)

I am trying to use samtools faidx to tkae this sequence but it doesn't work and keep returning me :

[fai_fetch] Warning - Reference 111-555 not found in FASTA file, returning empty sequence

After seeing some other users with the same problems, i tried to change the header name for it to contain no space (like ">scaffoldX") but it still doesn't work

Here is the exact command I type (my fasta file only contain the sequence of my scaffold) :

samtools faidx scaffold.fasta 111614-111868

Thanks a lot for your help.

Answer: Shouldn't it be "samtools faidx scaffold.fasta scaffoldX:111614-111868"?
Question: I have a file with around 20000 columns as gene names. I want to grep out the rpkm values for specific genes. Is there a way to grep out the column information?

> sample gene17 gene92 gene1 ... gene20000

> patient1 0.03569654 1.020565 0.0036522 ... 0.25247236

I only want gene72 for example, but it's not sorted in increasing order.

Thanks.

Answer: Hi,

in order to find the column number you can use:

head -n 1 file | tr '\t' '\n' | cat -n | grep gene72
With `head -n 1` , you get only the file's first line. With `tr` you replace the tab-separator by a new line. With `cat -n`, you print the input with line numbers on which you use finally `grep` to get the column of interest.

with the found number - let it be j - you can use cut:

cut -f 1,j file

Cheers,

Michael
Question: Hi everyone,

Sometimes I am asked to mention the most exciting moment or part of my research works, I often say that I was excited first time I know how to work with command lines in Linux. However I think might be there is a specific purpose behind asking this question. What would you reply if you asked so?

Thank you very much

Answer: As a part wet lab rat, part computer guy, the most exciting thing is actually the fact to be able to analyze the data you produced in the lab yourself, without being dependent on a second person (at least standard analysis). That directly leads to the most frustrating thing, which is realizing how overwhelmingly much I still cannot do, especially when it comes to advanced data analysis strategies and statistics.
Question: Hello!

Recentlly, I trimmed my poor tail .fastq dataset using Trimmomatic's SLIDINGWINDOW parameter, in which I chose threshold value 15~30(i.e SLIDINGWINDOW:4:15~SLIDINGWINDOW:4:30). after, I mapped them using STAR, but I found with the increase of threshold, totally, unique and multiple mapped reads were all decreased across my 5 test_samples.

So, I'm in doubt with this situation, should it's a normal sense that my totally mapped or unique mapped reads increase with the accompany of increased trimming threshold?

![Mapped results][1]

[1]: https://image.ibb.co/iHOBHk/trim_star_compare.png
I would really appreciate for your helping!

Answer: Sliding window trimming is not optimal. Because in this case there is an optimal algorithm (the Phred algorithm), I suggest you never use sliding window. The Phred algorithm is implemented in seqtk and BBDuk.

Different aligners handle low-quality bases differently. Quality-trimming will almost universally decrease the proportion of uniquely-mapped reads, from a given set of reads that map whether trimmed or untrimmed. Decreasing the rate of total mapped reads generally means that you are quality-trimming too aggressively, or simply using an aligner that is not sensitive enough. STAR is very fast, but not as sensitive as some other aligners; perhaps it does not handle the shorter reads resulting from quality-trimming very well, or perhaps it is very tolerant of low-quality bases (local aligners tend to be).

Anyway, the image would be more useful if you used optimal quality trimming and included thresholds below Q15, which is higher than is generally recommended for mapping (particularly for anything quantitative like RNA-seq). Q10 and Q5 are worth including, and from there you can bisect if you want the optimal value.
Question: Apologies if this has already been asked; I've looked around and haven't been able to get a clear, definitive answer to this question.

I have a number of interleaved paired-end FASTQ files that are sorted. I read in a (somewhat dated) GATK tutorial that the reads in these files should be sorted randomly (keeping the pairs together, one would assume) prior to alignment with BWA-mem, lest a bias be introduced. Is this correct?

Assuming it is, is there a tool that you know of / can recommend to sort interleaved paired-end FASTQs in such a way? I know there are scripts that can be written, but I'd prefer to stand on the shoulders of a (verified, preferably open-source) giant.

Many thanks!

Answer: `shuffle.sh` from [BBMap suite][1].

Description: Reorders reads randomly, keeping pairs together.

Usage: shuffle.sh in= out=

[1]: https://sourceforge.net/projects/bbmap/