cannin/biostars_qa huggingface dataset:
Question: Dear all,

I am struggling with what should probably be a trivial problem.
I have 2 seg files with the following structure:

chr start end feature.1 feature.2
1 1000 1200 0.5 0.001

chr start end feature.1 feature.2
1 1100 1400 0.7 0.004

Let's call them A and B. I want to obtain a unique file like this:

chr start end feature.1.a feature.1.b
1 1000 1100 0.5 NA
1 1100 1200 0.5 0.7
1 1200 1400 NA 0.7

If a region is covered by only A or B, then it would contain only the corresponding values and it would have NA in the other fields. Otherwise, I want to retrieve both of them.
This looks close to what bedtools -intersect does, and indeed I am trying to use the package Bedr to solve this task, but it seems to leave out regions covered only in one file.
Maybe GenomicRanges is another possibility.
If you have any hint on this it would be very helpful to me.
Thank you

Answer: I managed to deal with the issue using data.table::foverlaps. I report an example from another thread.

https://www.biostars.org/p/348933/#348990

The solution needs to be adapted but essentially solves the problem.
Question: Hi, I have 10 microarray data sets (each data set related to a disease) in which I already compiled to get differentially expressed genes (DEGs) in every individual data set. I want to extract common DEGs of all the datasets, is there any tool/R package/R function to do it.

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

Any help would be appreciated. Thanks!

Answer: If I understand it correctly, you simply want to find the unique gene identifiers (or probe ids) that are differentially expressed in all experiments? One simple way to do that in R would be (here for three sets):

a = c('gene1','gene3','gene5','gene7','gene9')

b = c('gene3','gene6','gene8','gene9','gene10')

c = c('gene2','gene3','gene4','gene5','gene7','gene9')

Reduce(intersect, list(a,b,c))

[1] "gene3" "gene9"

On a side note, if you are requiring a gene to be significantly differentially expressed in all experiments, that is a fairly tough threshold. Since experiments typically do not have the power to detect all genes that are truly differentially expressed, you are likely to miss some in each set due to randomness. You could, alternatively, require that a gene is significantly differentially expressed in some samples, and has a fold change in the same direction (or over some meaningful threshold) in others.
Question: I understand that unlike for ChIP-seq, with ATAC-seq MACS doesn't estimate the model, which means it doesn't determine distance 'd' between the tags and therefore it doesn't shift the tags by d/2 in 3´direction (and that is why with ATAC-seq the --nomodel --shift 0 parameters are set).

But if 'd' is not estimated then how can MACS slide a window of 2d across the genome to calculate peak enrichment?

Also, for ChIP-seq there is a control sample (input DNA or DNA pulled down with an unspecific antibody) that is used to calculate lambda local and determine peak enrichment. In ATAC-seq there isn't any control pulled down with unspecific antibody. What is that control in ATAC-seq?

Answer: MACS (for ChIP or ATAC) does not necessarily need a control. There is actually a section in the paper describing exactly the situation where no control is available. In this case, it determines the local lambda in certain windows. The paragraph starts with:

> Therefore, instead of using a uniform λBG, estimated from the whole genome, MACS uses a dynamic parameter, λlocal, defined for each candidate peak as (...)

In the absence of a control, naively one would estimate the background level of the experiment in a uniform fashion. That means if you have, say 30mio reads and you randomly threw them onto the genome, then every base (given a certain fragment length of the library) would have a coverage of x. Still, the genome coverage is never uniform (you will most impressively see once you analyze your first WGS sample, it is really a hilly landscape) due to differences in chromatin structure, PCR/GC bias, local copy number alterations etc. So instead of a genome-wide λBG, MACS checks the vicinity of the peak centers (up to 10kb) to estimate how prone this genomic region is to accumulate reads. In my understanding, as a typical ChIPseq experiment produces sharp peaks, the local environment should be depleted for enriched signals. Therefore, notable readcounts in the vicinity are an indication of a local bias. As a result, the peak enrichment needs to be penalized (down-corrected), as the region itself is prone to accumulate enrichment, irrespective of the protein target.
Question: sorry,

who know few words about mireap procedure, by googling I could not find even one syntax :(

Answer: Dear Angel, Hi

If you want some helps about script and syntax for running this software it is in the **README** file of the program and is as below, hope that helps:

Program: MIREAP (Reap miRNAs from deeply sequenced smRNA library)
Version: 0.2
Contact: Li Qibin
Bioinformatics department, Beijing Genomics Institute

**1. Introduction**
MIREAP combines small RNA position and depth with a model of
microRNA biogenesis to discover microRNAs from deeply sequenced
small RNA library.

**2. Installation**
You must have Vienna RNA Package (http://www.tbi.univie.ac.at/RNA)
installed on your computer and make sure that its perl interface
is accessible.

Copy mireap_0.2.tar.gz to a directory (/foo/bar) and unpack
it by command:
tar -zxvf mireap_0.2.tar.gz

Before running mireap, you need add path /foo/bar/mireap_0.1/lib to
environment variable PERL5LIB:
For csh/tcsh:
setenv PERL5LIB /foo/bar/mireap_0.2/lib
For sh/ksh/bash:
export PERL5LIB=/foo/bar/mireap_0.2/lib

**3. Usage**

mireap.pl -i -m -r -o

Options:

-i Small RNA library, fasta format, forced

-m Mapping file, tabular format, forced

-r Reference file, fasta format, forced

-o Directory where results produce (current directory)

-t Sample label (xxx)

-A Minimal miRNA sequence length (18)

-B Maximal miRNA sequence length (26)

-a Minimal miRNA reference sequence length (20)

-b Maximal miRNA reference sequence length (24)

-u Maximal copy number of miRNAs on reference (20)

-e Maximal free energy allowed for a miRNA precursor (-18)

-d Maximal space between miRNA and miRNA* (35)

-p Minimal base pairs of miRNA and miRNA*

-v Maximal bulge of miRNA and miRNA* (4)

-s Maximal asymmetry of miRNA/miRNA* duplex

-f Flank sequence length of miRNA precursor (10)

-h Help

Please convert your small RNA file into fasta format and append
sequencing frequence to sequence Id, just like this entry:
>t0000035 3234
GAATGGATAAGGATTAGCGATGATACA
(t0000035 is read_ID, 3234 is sequencing frequence)

The format of small RNA mapping file should be (delimited by tab or
space):
read_ID,chr_ID,start,end,strand(+/-)

You can make MIREAP run on the test data by execute comand:
perl ../bin/mireap.pl -i rna.fa -m map.txt -r ref.fa
.

.

**4. Output format**
MIREAP produce three files at each run.

*.gff
This file contains miRNA genes discovered by MIREAP, GFF3 format. For
GFF3 format, please refer to http://www.sequenceontology.org/gff3.shtml
Attribute 'Count' denotes the sequenceing frequence.

*.aln
This file contains sequence and structure of the pre-miRNA. Small RNAs
also are aligned to the precursor from which you can get more insights
into the maturation process of miRNAs.

*.log
This log file records parameters, start end time and other informations.
Question: Currently, I am trying to analyze some CAGE data from the FANTOM5 consortium (fantom.gsc.riken.jp and for the data http://fantom.gsc.riken.jp/5/data/).

Unfortunately, I am not that experienced in analyzing CAGE data. They do provide some processed data for both the TSS expression and enhancers and some processing software. However, I want reprocess these data myself and to start at least with the BAM files. So given CAGE BAM files, what would you suggest to use to get out both a sensible value of gene expression for each TSS (and Isoform) and the same for the enhancers?

For the former, usually I would simply take cuffdiff with some GENCODE annotation, but I am unsure whether this is just sensible for RNA-Seq and not for CAGE - or what would you use for it? And I still have no idea how to do it similarly with enhancers (under the assumption that I am having their annotation) - should I count the number of tags falling into these regions?. Also do you know how the "activity" of an enhancer is related to eRNA abundance.

Generally speaking, should I consider the mRNA abundance for a gene measured with CAGE in FANTOM5 as the level of transcription initiation or the the steady-state mRNA level after posttranscriptional regulation (by miRNAs)? Also are other RNAs like miRNAs or lncRNAs included?

It would be great if you can help me here.

Answer: CAGE detects promoters of coding and lncRNA, miRNA promoters, post-transcriptionally processed RNAs, eRNA and so on.

1. If you want to process CAGE BAM files yourself, here is the tool: http://www.ncbi.nlm.nih.gov/pubmed/25653163

2. CAGE detects post-transcriptionally processed RNA signals spread across exons, and they have very different characteristic features, here is the paper: http://www.ncbi.nlm.nih.gov/pubmed/24002785

3. CAGE detects miRNA promoters along with Drosha processing events, described here: http://www.ncbi.nlm.nih.gov/pubmed/26673698

In general, CAGE tags detected at the promoter is mRNA abundance, that is used to quantify gene expression level. CAGE tags at the promoters show nice correlation with RNA-seq. You can also read the Fantom5 papers itself.
Question: HI

Can some one give better explanation of primary and alternate genome assembly ?

The organism I'm working has only "primary assembly" can I use this assembly for my data analysis work?

Need suggestions!

Answer: NCBI's assembly model is described [on this page][1]. If a RefSeq assembly is available for your organism then you should use that for your (unspecified above) analysis.

[1]: https://www.ncbi.nlm.nih.gov/assembly/model/
Question: Hi All,

I am wondering if anyone estimated transcript-level counts from plate based scRNA datasets like smart-seq that capture full length transcripts instead of 3' ends. I tried to use salmon but the estimated counts/TPMs seems to be highly inflated.

For example, for this transcript, the output from a single cell looks like this.

Name Length EffectiveLength TPM NumReads
ENST00000344843.7 721 250 247.016 471.945

When I looked the same cell bam file in IGV, that transcript has only around 50 reads mapped. Salmon `NumReads` is 471. This happens for a lot of transcripts.

![enter image description here][1]

I am wondering why the values are inflated ? One potential reason could be due to the default scaling factor used, as for bulk-rna where the total counts tend to be in millions.

I would like to know if I anyone estimated the transcript counts from scRNA before. There are tools like [Alevin][2] but they seem to work only with 3' enriched droplet based methods.

Thanks,

Goutham A

[1]: https://i.ibb.co/202n9qQ/salmon.png
[2]: https://salmon.readthedocs.io/en/latest/alevin.html

Answer: Since this is full length without UMIs, salmon would be the tool to use here. Out of curiosity, how did you obtain the pile-up? How do you know those are all reads coming from this transcript? Is it a single transcript gene? You can always align the reads with STAR and get a bam file using `—quantMode TranscriptomeSAM` if you want quantification from those alignments.
Question: Dear All

I have a problem that we can not solved for a while in installing bioawk. I follow the general pipeline to install bioawk on ubuntu server:

- clone from git hub to the Download directory
- make

After that it is where we are stuck:

When I run the command:

mv awk book

I got error awk does not exist! When I install gawk and try `mv gawk` from `/usr/bin` to `home/arsalan/download/bioawk`, it works but when I `cp bioawk` to `usr/local/bin` and try to run it if fails and since gawk has already moved from `/usr/bin`. By typing `gawk` I got the error it is not installed! and typing `bioawk` I got error command not found!

Many thanks in advance

Answer: bioawk now compiles to the executable called `bioawk` and not `awk`

it used to compile to awk and some instructions may still refer to that and instruct you to rename it
Question: I am wondering what the best practices are analyzing host-parasite mixed transcriptomes. If I have mRNA-Seq data from infected host tissues that contain both host and parasite genomes, which would be the best way to tease them apart? Mapping to a combined reference (host + parasite), or mapping to the host genome first, collecting unmapped reads and mapping to parasite genome? What would be the implications of these options? When I get the real data, I could run tests with both approaches, but maybe someone already have some thoughts on this here?

Answer: You should use `bbsplit.sh` from BBMap suite for this purpose. It will allow you to bin the reads to multiple genomes (and handle multi-mappers in intelligent ways). See this [thread][1].

[1]: http://seqanswers.com/forums/showthread.php?t=41288
Question: Hello all,

I am looking to test an alignment-free phylogenetic tree building algorithm I wrote. It can perform both gene and species trees. I have already tested it on a single gene primate tree, but I need some more data to further characterize the algorithm. I know there is a lot of data on TreeBASE, but I am having a hard time pulling data down. Additionally, I am generally unaware of which trees are considered well-resolved.

Any info would help greatly

Answer: You might want to use data sets already used in other papers on alignment-free comparisons. [Here][1] you can download the data from andi (shameless self-plug). I also have the roseobacter data set from the [spaced words][2] paper. Send me a mail, if you are interested.

[1]: https://github.com/EvolBioInf/andi#data-sets "andi data"
[2]: http://spaced.gobics.de/ "spaced words"
Question: I am having trouble using the crossmeta package for analyzing GEO datasets. The library downloaded fine and I got through the first few steps without any issues:

library(crossmeta)

# specify where data will be downloaded
data_dir <- file.path(getwd(), "data", "LY")

# gather all GSEs
gse_names <- c("GSE9601", "GSE15069", "GSE50841", "GSE34817", "GSE29689")

# gather Illumina GSEs (see 'Checking Raw Illumina Data')
illum_names <- c("GSE50841", "GSE34817", "GSE29689")

# download raw data
# get_raw(gse_names, data_dir)

But once I tried the next command:

# this is why we gathered Illumina GSEs
open_raw_illum(illum_names, data_dir)

I got this error:

Error in open_raw_illum(illum_names, data_dir) :
could not find function "open_raw_illum"

Can anyone help me figure out how to get it to recognize the command? Here is the link to the guide I am trying to follow: https://bioconductor.org/packages/devel/bioc/vignettes/crossmeta/inst/doc/crossmeta-vignette.html

Answer: Okay, thank you. I did look through all of the commands for the package, and that one is not listed, so I have no idea where it came from. I'll reach out to the author and see if I can get any insight.

- [GitHub issue 8](https://github.com/alexvpickering/crossmeta/issues/8) Code "open_raw_illum" not working
Question: I have downloaded the Level_3 data for CNV (SNP Array) for Breast Cancer from the Bulk Download section of TCGA.

The file names look like this:

MOHEL_p_TCGA_271_274_275_N_GenomeWideSNP_6_A01_1320320.hg18.seg.txt

The data in the "Sample" column within the file is the same.

I was expecting a barcode like `TCGA-A2-A0D0-01`.

How do I match the data to the patient?

Thanks for helping

Stephanie

Answer: I think you should download the data using the Data Matrix. When you download the data, you automatically download files like `FILE_SAMPLE_MAP.txt` & `file_manifest.txt`. You will get the barcode along with the filename in these files.
Question: I would like to download the mRNA expression z-scores for all genes (in all patients, with all cancers) from CBioPortal. I've found the download page, but I have to enter a list of genes I want the z-scores for - I don't want particular genes, I want to download all of them.

This data is normally available from the TCGA, but in this case I understand the z-scores are calculated from TCGA data by CBioPortal, and I'd rather not recalculate them myself. Is it possible to do a bulk download of all z-scores from CBioPortal?

Thanks for any help,

Stephanie

Answer: cBioPortal, unfortunately, does not support such bulk-download requests. As you also mentioned, you can download these data sets from TCGA; since you want to work with higher-level data, I suggest you download the data from Firehose analysis results instead of using TCGA's raw data download site:

http://gdac.broadinstitute.org/

More information is available on this web site, but you might want to check out the firehose-get tool, which will let you script your downloads in an easy manner.
Question: Hi all.I am dealing with some miRNA data sequenced by illumina Htseq 2000.And I confuse the result of the mirdeep2 software.
First is the result of the quantifier.pl.In the miRNAs_expressed_all_samples_time.csv has entries show below.

hsa-miR-1 549.00 hsa-mir-1-1 549.00 420.00 129.00 19.07 23.63
hsa-miR-1 564.00 hsa-mir-1-2 564.00 432.00 132.00 19.62 24.18

The same miRNA with different precursor.I am not sure how to deal with these when I try to use the R package DEGseq to detect
DE miRNA.Treat them as two different entries or just plus them if these entries have no intersection.

Second,the result of less output.mrd |grep '>'|wc -l is 1789.I think it hsa 1789 precursors,and the result.csv file just have 785 novel and known miRNAs, I don't know how to explain the different between these two numbers.
Your reply will be appreciated.
Cheers.

Answer: Different miRNA precursors can produce the same mature form.

In such cases, it depends on the type of your downstream analysis.

If you are interested in precursor expression then these two (and all others that fall in this category) should be considered as separate features.

If your analysis focuses on mature miRNAs (for target identification, etc) then these two matures should be considered as one, since they have the same mature sequence (same targets, etc).

Depending on your miRDeep 2 settings (if you divide multimaps or not), you can just add the expression of those two mature miRNAs into one, or keep one with just their average.

Cheers,

IV
Question: I have a bunch of viral genbank assembly identifiers such as this:
`['GCA_001041635.1', 'GCA_002958635.1', 'GCA_000915755.1', 'GCA_001041575.1']`

How can I get the actual taxonomic names for these?

My process right now is cumbersome:
1. Look up assembly ID on NCBI
https://www.ncbi.nlm.nih.gov/assembly/GCF_001041635.1/

2. Click on Related Information -> Taxonomy
https://www.ncbi.nlm.nih.gov/taxonomy?LinkName=assembly_taxonomy&from_uid=4780598

3. Click on the taxon under "Links from Assembly":
https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1500814

4. Copy "Lineage":
Viruses; Duplodnaviria; Heunggongvirae; Uroviricota; Caudoviricetes; Caudovirales; Siphoviridae; Pahexavirus; Propionibacterium virus PHL095N00

**How can I automate this process?** I made a crude webscraping function but I get locked out and it's not reliable nor is it robust.

Essentially I'm looking for one of the following approaches:

1. Python module that is either pip or conda installable that has a function: `f(genbank_assembly_id) -> taxonomic_lineage_string` which is probably unlikely.

2. A flat file that has at least 2 columns where one is Genbank Accession ID and another is lineage.

Answer: Using [EntrezDirect][1]:

$ esearch -db assembly -query "GCA_001041635" | elink -target taxonomy | efetch -format native -mode xml | grep ScientificName | awk -F ">|<" 'BEGIN{ORS=", ";}{print $3;}'
Propionibacterium phage PHL095N00, Viruses, Duplodnaviria, Heunggongvirae, Uroviricota, Caudoviricetes, Caudovirales, Siphoviridae, Pahexavirus, Propionibacterium virus PHL095N00,

I will leave crafting a loop to get output described in #2 up to you.

[1]: https://bit.ly/entrez-direct
Question: I have to calculate the Shannon entropy for a given list of sequences in a fasta file. Recently I came across a program called rnaentropy which solves my problem but the issue is I cannot run sequences in batch. In the documentation the input is given as a sequence of strings. The parameters in the documentation is as follows:

./RNAentropy "sequence" -s sequence -t temperature -e energyModel(default 2004) -d delta_Temp (compute structural entropy using <E> = RT² * d/dT ln(Z(T)) ) -c (centered) -z energy_is_zero [0|1] (default 0) -v (verbose, extended output)
./RNAentropy -h (detailed help)

Input paramters are:
<sequence> : If FIRST argument is a valid nucleotide sequence it will be used input
-s <sequence> : Alternative flag for input sequence
-t <temperature> : Temperature in ºC (default is 17ºC)
-e <energy_model> : Thermodynamic energy model used. Valid values for energy model are 1999, 2004 and 2007 (resp. Turner'99, Turner'04 and Andronescu '07) (default is 2004)
-d <delta_T> : Temperature variation used for estimating d/dT ln(Z(T))
NOTE: If this parameter is provided, H is computed estimating <E> = RT² * d/dT ln(Z(T))
-c : (Use only in combination of -d) Use the centered version for estimating <E> = RT² * d/dT ln(Z(T))
-v : Output includes method for computing H and the name of the ouput parameters
-z <energy_is_zero>: [0|1] If value is 1, energies are set to 0, ouput is structural entropy for the uniform case

Can anyone please suggest me a way to do it for a multi-fasta? Is there a way to iterate in a shell script?

Answer: With GNU Parallel (https://www.biostars.org/p/63816/) you would do:

cat in.fasta |
parallel --pipe -N1 --recstart '>' --rrs 'read a; echo Doing "$a"; ./RNAentropy -s $(tr -d "\n")'

From [GNU Parallel Tutorial: Informational](https://www.gnu.org/software/parallel/parallel_tutorial.html#Informational):
> `--recend str` Record end separator for `--pipe`.
> `--recstart str` Record start separator for `--pipe`.
Question: I have 7000 genes and their proteins as well as the genome of a bacteria am working on, I want to convert these files into either genebank format or embl, my problem is I don't have any scripting skills, I tried using an online tool (http://genome.nci.nih.gov/tools/reformat.html) but its out put appears to lack some information.

Does any one know a tool or a way I can convert these sequences?

Answer: No need to develop tools to do this. Many are publicly available for such a common task.

The standard for many years has been Emboss' SeqRet. An online version is here, but I would consider installing the suite if this is something you need to do often. The command line version is as simple as seqret <in> <out>. Wrap it in a loop.

BioPython's SeqIO module can also do this, albeit with a bit more (basic) programming. I'm sure there are equivalents in BioPerl, BioRuby, and via Bioconductor for R.

This kind of task is day one bioinformatics, and the skills required are easy to learn and very straightforward. It's as simple as navigating to a folder and running a program, possibly within the simplest of loops depending on how your data is organized. You are already converting 7000 sequences; it would make sense that you learn about the plethora of resources available to you developed over the past couple decades. You'll also save a ton of time in the future!