cannin/biostars_qa huggingface dataset:
Question: I have downloaded the miRNA Isoform expression quantification data (`mirbase21.isoforms.quantification.txt`) from TCGA. The data looks like below:

miRNA_ID	read_count	miRNA_region
hsa-let-7b	2	precursor
hsa-let-7b	1	precursor
hsa-let-7b	1 precursor
hsa-let-7b	58	MIMAT0000063
hsa-let-7b	173	MIMAT0000063
hsa-let-7b	5723	MIMAT0000063
hsa-let-7b	26947	MIMAT0000063
hsa-let-7b	1	stemloop
hsa-let-7b	1	MIMAT0004482
hsa-let-7b	2	MIMAT0004482
hsa-let-7b	129	MIMAT0004482
hsa-let-7b	401	MIMAT0004482

And based on `miRBaseConverter` and information from `miRTarbase` I have the mature regions 3p and 5p information with Accession like below:

miRNAName_v21	Accession
hsa-let-7b-5p	MIMAT0000063
hsa-let-7b-3p	MIMAT0004482

So, based on above information I can sum the counts of mature 3p and 5p. And the main precursor also I can do, but what are stemloops? Do I need to include that with precursor or need to exclude that from analysis?

thanq

Answer: The the maturing of a miRNA is a multi-step process. The extra material at the 3' and 5' end is removed first, and then the loop, which acording to the TCGA pipeline description defined for the purposes of this pipeline as:

> 4. stemloop, from 1 to 6 bases outside the mature strand, between the mature and star strand

In the TCGA pipeline, reads are classified heirarchically.
Priority is given to reads mapping to mature miRNA, so if it aligns to something else as well as mature, it will be counted as mature (that actaully wouldn't be the way I would do, but never mind). Then to reads mapping to the pre-cursor, and then to the stemloop.

This explains the high levels of counts in the mature, but not the pre-cursor, despite the mature being contained within the precursor.

As you are probably intersted in the biologically active molecule, I would ignore everything other than the mature.
Question: Hi,

I want to use fastQC for quality control checks on nanopore sequencing reads. When I ran fastQC it ran a long time without finishing or given an error. I stopped the job after several hours and got an error: "java.lang.OutOfMemoryError: Java heap space"
What is the reason for it? Is fastQC good for nanopre reads? (I read that fastQC added support for processing files of Oxford Nanopore reads), Which other tools are recommended for quality control checks on nanopore sequencing reads?

Thanks

Answer: Here are a few tools I've used..I think the question comes down to what exactly do you mean by QC (what stats are you trying to find out):

[minion_qc][1]

[IONiser][2]

[Nanoplot..as seen here on biostars][3]

[1]: https://github.com/roblanf/minion_qc
[2]: http://bioconductor.org/packages/devel/bioc/html/IONiseR.html
[3]: https://www.biostars.org/p/256345/
Question: Hi,

Anyone familiar with Jbrowse genome browser here? I wonder, can Jbrowse be linked to other party software such as Galaxy or Taverna? I found a tool named ABrowse (http://www.abrowse.org/), and seems like ABrowse can do that. Anyone?

Answer: I haven't seen ABrowse in action, but I looked at what are some features that ABrowse claims to have and here are some ideas about how JBrowse might offer similar functions

- On-the-fly submit selected entries to Galaxy and WebLab for further analysis.
You can add custom javascript callbacks to JBrowse to allow it to perform a variety of actions and so you could possibly use this to submit jobs to Galaxy.
You might ask the JBrowse mailing list to see if they have examples of Galaxy integration.
- Compatible with BioMart engine for bulk data retrieval.
I know that BioMart has a SPARQL query interface, and since JBrowse supports SPARQL backend, you could possibly use JBrowse to retrieve and display data from BioMart.
- SOAP-based web service interface, fully compatible with Taverna.
JBrowse has no backend server side code, it's a 100% client based javascript app, so it doesn't really offer web-services like ABrowse might. There are some other software systems such as Intermine that have better integration with JBrowse and web services though.
- Easy to install, configure and customize for different user requirements.
I think JBrowse is also easy to install, configure, and customize (extensive customizations possible via javascript! no complex server side code :)

Feel free to ask the gmod-ajax mailing list for jbrowse topics, beginner or advanced!
Question: Hello,

does anybody of you has ever tried to generate null model of DNA sequences in mouse using the gkmSVM package? It works perfectly for human but for mouse it does not. I am wondering if any of you have ever used this package for mice and had the same problem.

The error problem when I use the genNullseqs function is:

library(gkmSVM)
library(BSgenome.Mmusculus.UCSC.mm10)
library(BSgenome.Mmusculus.UCSC.mm10.masked)
library(IRanges)

genome=BSgenome.Mmusculus.UCSC.mm10.masked

fileBedBreaks="Rep1.intersec.Rep2.cov.major2.sort.uniq.bed"
fileFastaPos="Rep1.intersec.Rep2.cov.major2.all.sort.uniq.bed.pos.fa"
fileBedNeg="Rep1.intersec.Rep2.cov.major2.all.sort.uniq.Random.gkmSVM.bed"
fileFastaNeg="Rep1.intersec.Rep2.cov.major2.all.sort.uniq.Random.gkmSVM.fa"
genNullSeqs(inputBedFN=fileBedBreaks,nMaxTrials=5,xfold=2,genome=genome,
outputPosFastaFN=fileFastaPos,outputBedFN=fileBedNeg,outputNegFastaFN=fileFastaNeg)

`Error in normalizeDoubleBracketSubscript(i, x, exact = exact) : subscript "TRF" matches no name`

TRF is related with the tandem repeats because it means Tandem Repeats Finder. But what does it means the problem??

Thanks in advance for any reply

Answer: I hit this same problem and figured out that it's because the BSgenome.Mmusculus.UCSC.mm10.masked on Bioconductor (as of May 2020) only has two masks included (AGAPS and AMB), whereas previous genomes included the TRF mask as well, which is what gkmSVM is looking for.

Thankfully, on GitHub, Zhang (Frank) Cheng has added the TRF mask to the mm10 genome at [https://github.com/biomystery/BSgenome.Mmusculus.UCSC.mm10.masked][1]

I had to uninstall the Bioconductor version of mm10 and reinstalled using `remotes::install_github("biomystery/BSgenome.Mmusculus.UCSC.mm10.masked")`

[1]: https://github.com/biomystery/BSgenome.Mmusculus.UCSC.mm10.masked
Question: I am using Eage to phase plink input. I maintain --keep-allele-order throughout my workflow. Eagle outputs .haps (Oxford phased haplotype file), which I convert to VCF with shapeit. I did some sanity checks and find that this process reverses REF/ALT designation and genotypes for every single one of my snps and indels. Comparing alt allele frequency to reference strongly suggests that I do not have a build issue, so the workaround was initially to do

`bcftools norm --check-ref -s -f hg19.fa input.vcf > output.vcf`

to flip everything back after converting to vcf with shapeit, which perfectly flips every single one of my snps, but does not seem to do the same for indels. Is there a way to do this for indels as well. My alternative thought is just to use bcf as input for eagle, which should avoid the Plink/Oxford formats completely.

Answer: You should be using plink 2.0 instead of 1.9 as much as possible whenever allele order matters. Otherwise, if you forget --keep-allele-order even once, you're in trouble; and there are some plink 1.9 commands (like --linear/--logistic with interaction terms) which need to NOT be run with --keep-allele-order.

You can use plink2 [--ref-allele][1] to reset the REF alleles in the manner you want; see in particular the second bullet point in the documentation.

And as for *every* variant having REF/ALT reversed, that sounds like you did not import the initial file correctly. Recent plink 2.0 builds make that particular mistake a lot less likely by forcing you to declare whether an Oxford-format file has REF first or last.

[1]: https://www.cog-genomics.org/plink/2.0/data#ref_allele
Question: My RNA-seq read count table consists of samples fro two tissues (root and shoot), different time-points (1hr, 3hr, 12hr, 1day etc). and stress and no-stress conditions all in one table.
For example:

`root_ctrl_rep1 root_ctrl_rep2 root_1hr_rep1 root_1hr_rep2 root_3hr_rep1 root_3hr_rep2 shoot_ctrl_rep1 shoot_ctrl_rep2 shoot_1hr_rep1 shoot_1hr_rep2 shoot_3hr_rep1 shoot_3hr_rep2 `

As of now, I am interested only in a basic DEG analysis between the following pairs

`root_ctrl vs root_1hr`

`root_ctrl vs root_3hr`

`shoot_ ctrl vs shoot_1hr`

`shoot_ctrl vs shoot_3hr`

----and so on

How do go about this ?.

> library("DESeq2")

>count_table<-read.table("counts.txt", header=T, row.names=1, sep="\t")

>exp_design<-data.frame(row.names=colnames(count_table), condition = c("root_ctrl","root_1hr","root_6hr","root_12hr","shoot_ctrl","shoot_1hr","shoot_6hr","shoot_12hr","shoot_day1","root_ctrl","root_3hr","root_6hr","root_12hr","root_day1","shoot_ctrl","shoot_1hr","shoot_3hr","shoot_6hr","shoot_12hr","shoot_12hr","shoot_day1","root_1hr","root_3hr","root_day1","shoot_3hr","shoot_3hr"))

>conditions=exp_design$condition

> ds <- DESeqDataSetFromMatrix(countData=count_table, exp_design, formula(~condition))

I am not sure if the above is right and what to do next. Any suggestion is greatly appreciated. Thanks

Answer: You'll probably want `~0+condition`, since I'm not sure if DESeq2's extended model matrix will be used in your model or not (I presume not).

Then for each comparison: `res = results(ds, contrast=c("condition", "root_1hr", "root_ctrl"))`. Note that you'll want to include `lfcShrink()`, since that's not automatically done anymore.
Question: I'm trying to get some `entrez` ids from some gene names using [biomaRt][1]

In my example below I have 2 genes `Igha` and `Mlc1` for mus musculus

My version of biomaRt is `biomaRt_2.38.0`

library("biomaRt")
mart <- useMart(biomart = "ensembl", dataset = "mmusculus_gene_ensembl")
getBM(attributes=c("ensembl_gene_id","mgi_symbol","entrezgene"), values = c("Igha", "Mlc1"), bmHeader = T, filters = "mgi_symbol", mart = mart)
# Gene stable ID MGI symbol NCBI gene ID
#1 ENSMUSG00000095079 Igha NA
#2 ENSMUSG00000035805 Mlc1 170790

I know that some Ensembl gene ID do not have Entrez gene ID https://www.biostars.org/p/275619/

But in this case, I'm able to find the Entrez gene ID from Ensembl

http://www.ensembl.org/Mus_musculus/Gene/Summary?db=core;g=ENSMUSG00000095079;r=12:113254830-113260236

Which is `Refseq Gene ID 238447`

Maybe https://www.biostars.org/u/7019/ will have an answer :)

[1]: https://bioconductor.org/packages/release/bioc/html/biomaRt.html

Answer: Using https://www.biostars.org/u/18713/ 's comment

Here is my trick using [esearch][1]

biomart_genes <- getBM(attributes = c("external_gene_name", "entrezgene"), filters = "mgi_symbol", values = gene_list, bmHeader = T, mart = mart)
colnames(biomart_genes) <- c("hgnc_symbol","entrez_gene_id")

for (i in 1:nrow(biomart_genes)){
if (is.na(biomart_genes[i,"entrez_gene_id"])){
term <- paste0(biomart_genes[i,"hgnc_symbol"]," [gene] and 10090 [taxID]")
biomart_genes[i,"entrez_gene_id"] <- uid(esearch(term, db = "gene", rettype = "uilist", retmode = "xml", retmax = 1))
}
}

[1]: https://www.rdocumentation.org/packages/reutils/versions/0.2.2/topics/esearch
Question: I am doing gene ontology enrichment with cluterprofiler package, and I got the result but when I plotting barplot (10000 descriptions) it's difficult to read or identify the name in the plot, is there any method to remove those descriptions which don't have any gene name.

Command Line

>ggo <- groupGO(gene = UPT2,
+ OrgDb = org.Hs.eg.db,
+ ont = "MF",
+ level = 3,
+ readable = TRUE)

>barplot(ggo, drop=TRUE, font.size = 10,showCategory=10000, title="Up Treated2 Biological Process"

Answer: > is there any method to remove those descriptions which don't have any gene name.

```
barplot(ggo, drop=TRUE, ...)
```

You already know how to do it using `drop=TRUE` parameter.

> I plotting barplot (10000 descriptions)

There is no 10000 terms in level 3 of MF.

> it's difficult to read or identify the name in the plot

export the plot to file by setting large value of width and height should works.

```
p <- barplot(...)
ggplot2::ggsave(p, filename = your_file, width=set_a_large_value, height=set_a_large_value)
```
Question: Hi there!

Is it possible to use categorical traits (such as ethnicity) in WGCNA analysis? If so, shall we encode them with 0s and 1s? what if the trait is not binary and has several levels? Perhaps use 1, 2, 3,4?

Thank you in advance

Answer: You can use categorical variables, if it is binary, then you can expressed it as 0 and 1. If it has more levels then you can separate each level as an individual binary variable.

[This link may help][1]

[1]: https://peterlangfelder.com/2018/11/25/working-with-categorical-variables/
Question: I am using this as a learning exercise, I realize there are existing tools to determine the end result. The goal is to gain experience with Python and nuances of syntax.

Given that I have a fasta file I am parsing, and want to find non-overlapping ORFs in all possible frames (0,1,2). My code I have written thus far is below:

def getORF(seq, treshold, start, stop):
start = ["ATG"]
stop = ["TAA","TAG","TGA"]
start_codon_index = 0
end_codon_index = 0
start_codon_found = False
for frame in range(3):
q=frame
for i in range(q, len(seq), 3):
current_codon = seq[i:i+3]
if current_codon in start and not start_codon_found:
start_codon_found = True
start_codon_index = i
e=start_codon_index+1
if current_codon in stop and start_codon_found:
end_codon_index = i
d=end_codon_index+3
length = end_codon_index - start_codon_index + 3
if length >= treshold * 3: #treshold is protein length but length is nuc
b=length//3 + 3
start_codon_found = False
print("aa %i, nuc %i, frame %i, coord %i:%i\n" % (b, length, q+1, e, d))

However, when I run the following (myfile contains a single sequence record, but goal is to use multiple fasta record) it only returns results for the first 2 frames.

for record in SeqIO.parse(myfile,'fasta',IUPAC.unambiguous_dna):
seq=record.seq
name=record.id
print("%s\n" % name)
getORF(seq,2,start,stop)

I know there are ORFs in frame 3, and I can find start and stop codons using the following scripts:

def test(seq, start, stop):
start = ["ATG"]
stop = ["TAA","TAG","TGA"]
start_codon_index = 0
end_codon_index = 0
q=2
for i in range(q, len(seq), 3):
current_codon = seq[i:i+3]
if current_codon in start:
start_codon_index = i
print("Start at %i " % start_codon_index)

def test2(seq, start, stop):
start = ["ATG"]
stop = ["TAA","TAG","TGA"]
start_codon_index = 0
end_codon_index = 0
q=2
for i in range(q, len(seq), 3):
current_codon = seq[i:i+3]
if current_codon in stop:
stop_codon_index = i
print("Stop at %i " % stop_codon_index)

for record in SeqIO.parse(myfile,'fasta',IUPAC.unambiguous_dna):
seq=record.seq
name=record.id
print("%s\n" % name)
test(seq,start,stop)

for record in SeqIO.parse(myfile,'fasta',IUPAC.unambiguous_dna):
seq=record.seq
name=record.id
print("%s\n" % name)
test2(seq,start,stop)

I can't figure out why I am getting no results for frame 2. I have even tried throwing in 4 as my range for frame, and instead I get a result for a non existent frame 4 and no results for frame 3.... Clearly I am missing something obvious! Thanks in advance for any help!

Answer: You have to reset your `start_codon_found` variable to `False` in each `frame` loop. Otherwise it can be that you find a start codon in frame 2 but no more stop codon. Then, in the next `frame` loop, your `start_codon_found` varible is still true. Try this:

...
for frame in range(3):
start_codon_index = 0
end_codon_index = 0
start_codon_found = False
for i in range(frame, len(seq), 3):
....

This makes it working for me.
Question: Hi everyone,

I have RNA-Seq data from 22 samples and 3 batches. I used bioconductor to generate the RPKM values. I want to see if there is any batch effect in my dataset. How do i use RPKM matrix as an input to perform PCA ? Any help is appreciated.

Thanks,

apoorva

Answer: I usually make my PCA plots with VST data from DESeq2 in Bioconductor. The following might work for RPKM... but I haven't tried it.

Put your VST matrix (try RPKM matrix) instead of “MATRIX” below (you might have to transpose it using the "t()" function):

pca= prcomp( MATRIX , center=T, scale=T)

pdf("pca.pdf")
plot(pca$rotation[,1],pca$rotation[,2], xlab = "PC1", ylab = "PC2")
text(pca$rotation[,1],pca$rotation[,2], row.names(pca$rotation), cex=0.5, pos=4)
dev.off()
Question: I am pretty new to bioinformatics, and just installed bash and ubuntu on my PC. When I use the 'ls' command for my home directory, it shows only a few of the files and folders. It also neglects to show my genome vcf file. Terminal also does not recognize that there is a vcf file present when i try to read it in using vcftools....I'm thinking these problems are related. I've had success on the Mac computers in my lab on campus. Please help!

Thanks.

Answer: I figured it out, thanks everybody! Some of my files didn't have permission. So I used the ls -a to find out which files those were, then used sudo chmod 777 filename to enable access to them via terminal.
Question: I have used `featureCounts` to get the read counts one a gene level with the following command:

featureCounts -a gencode.gtf -F GTF -p -s 2 -T 8 -o sample.txt sample.sorted.bam

I'm interested in extracting the counts on a transcript level. My data is paired end and reverse stranded. I'm a bit confused with command I have to use to get the counts on transcript level [featureCounts manual][1]

Can anyone please tell me how to get the counts on transcript level? thanq

[1]: http://bioinf.wehi.edu.au/featureCounts/

Answer: I am one of the featureCounts authors.

Do you want do conduct analyses for genes, or do you need individual results for all the isoform variations of each gene? The "isoforms" tend to be called "transcripts" in the RNA-seq literature. The featureCounts help page you link to is designed to generate gene-level counts. We have never recommended the use of featureCounts to produce transcript-level counts.

A gene level analysis is easiest to interpret and is appropriate for most analyses where you want to interpret the results in terms of pathways, GO terms, KEGG and so on. featureCounts will give you gene level counts very quickly and easily. Alternatively, featureCounts can give you counts for each exon in each gene.

If you do really want to quantify all the different transcripts (isoforms) for each gene, then this cannot be done by featureCounts or by any other read counting software. While featureCounts can be run with transcript annotation, the results won't be particularly useful. Since the transcripts are heavily overlapping, with many transcripts sharing the same exons, featureCounts cannot tell which transcript to assign each read to. It therefore has to either discard all the multi-assigned reads or alternatively assign them multiple times to all the overlapping transcripts -- and either choice is bad.

If you really do want the (considerable) extra complexity of transcript-level analysis, then I would run Salmon and read the results into edgeR using `catchSalmon`. Salmon can assign reads to transcripts on a probabilistic basis. Or alternatively, do an exon-level analysis. That allows you see what is going on, although the results are not assembled into transcripts. Exon-level may be a safer approach if your data has expressed transcripts that have not been previously annotated.
Question: I have a series of BWA index files ( *.amb *.ann *.bwt, etc) but the original reference fasta sequences are not available. Assuming that it's a complete burrows-wheeler index, the complete reference should be encoded in those indices.

Does anyone know if/how the fasta can be re-extracted from them?

Answer: This is an interesting question. It definitely should be possible, so I guess the question may be how to do it with existing tools and the least amount of pain.

I looked into it, and it seems that from the [original bwtsw code][1] (that bwa relies on, particularly for the indexing steps), the "packed" fasta file (`pac`) is the key file to use. This is supported by looking at the order of the indexing substeps of `fa2pac`, `pac2bwt`, and finally `bwt2sa`.

So as an informed guess I searched for "pac2fasta", and (surprisingly to me) found an existing utility in the [TMAP package][2]. I believe, based on some quick tests, that the binary format for their `pac` file is the same as bwa uses (and bwtsw), aside from some embedded version numbers. The only catch is that TMAP has a binary annotation file (`$ref.tmap.anno`) that stores the chromosome names and lengths, whereas bwa uses a plaintext annotation file (`$ref.ann`). So if this is really what you want to do, it looks like you'll have to hardcode things (that's what I did to test with a small single chromosome reference) or write some code to make a TMAP annotation file from a bam header or bwa `.ann` file (by working through `tmap_refseq_write_header` in [tmap_refseq.c][3]).

[1]: https://github.com/mruffalo/bwt-sw
[2]: https://github.com/iontorrent/TMAP/
[3]: https://github.com/iontorrent/TMAP/blob/master/src/index/tmap_refseq.c
Question: Hello,

Any one who can tell me which software can draw this figure?

Thanks in advance!

Yue

Answer: It's called a lollipop plot. A lot of packages exist that can draw plots like this, and what you pick depends on what you need. For example, https://github.com/freezecoder/mutsneedle is a nice R package, but there are things in there that you cannot change and that is extremely irritating. cBioPortal has one too, but it's a tool where you will need to upload your data and generate the plot manually.
Question: Hi all,

I am looking for an online repository for BAM files for teaching and testing purposes. It doesn't matter the genome they are mapped onto, but ideally they should be no bigger than 100MB. ENA (the European Nucleotide Archive) has some projects containing BAM files, but I cannot find them since this is not a queriable feature (the standard raw file is obviously FASTQ).

Can you help me? I promise you will be credited in all my lessons forever!

PS I know this post is a bit of a repetition of a previous one from 2014 ([link][1]), but the links shown there are either expired or contain files that cannot be downloaded (try it).

[1]: https://www.biostars.org/p/45320/

Answer: In the 1000 Genomes there are plenty of BAM file, check the http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/ location you can subselect alignments to chr11, chr22, or exomes.
Question: I have a gbk of a draft bacterial genomes (meaning multiple contigs).
I want to extract some info about specific genes, and I get an error I cannot understand.

Parse the gbk

record_dict = SeqIO.to_dict(SeqIO.parse("SpeciesA.gbk", "genbank"))
My file looks:

for i in record_dict:
... for f in record_dict[i].features:
... print f.qualifiers
{'locus_tag': ['AA_03640'], 'gene': ['ftsH_4']}
{'locus_tag': ['AA_03640'], 'inference': ['ab initio prediction:Prodigal:2.6', 'similar to AA sequence:UniProtKB:P37476'], 'codon_start': ['1'], 'EC_number': ['3.4.24.-'], 'transl_table': ['11'], 'product': ['ATP-dependent zinc metalloprotease FtsH'], ........

Now if I try to get out the contig name where the gene AA_03640 is in, as follow:

for i in record_dict:
... for f in record_dict[i].features:
... if f.qualifiers['locus_tag'] == 'AA_03640':
print(i)
But I get the following error:

Traceback (most recent call last):
File "", line 3, in
KeyError: 'locus_tag'

Any help to understand where I am wrong?
Thanks

Answer: Not all feature types will have a locus tag - you probably only care about CDS or gene features. Even then, not genes will have a locus tag (although generally within a record they either all do, or all don't).

Finally if there is a locus tag, then ``f.qualifiers['locus_tag']`` will return a list (since in general the annotation qualifiers can appear more than once), so you'd want to check using ``'AA_03640' in f.qualifiers['locus_tag']`` or similar, e.g.

for i in record_dict:
for f in record_dict[i].features:
if f.type == "CDS":
if 'AA_03640' in f.qualifiers.get('locus_tag', []):
print(i)

See also http://www.warwick.ac.uk/go/peter_cock/python/genbank which does something similar to what I think you want to achieve.