cannin/biostars_qa huggingface dataset:
Question: Aloha,

I am having trouble figuring out how to remove everything after the last '_' in the sequence headers of a fasta file.

I would like this following series of headers

>ART01B_100_M7_ID100005_1
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325_189
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005_46
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

to look like this:

>ART01B_100_M7_ID100005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

because there are various '_' in the sequence headers, this command, for example, won't work:

cat Alt_MACSE_Output.fasta | awk -F _ '/^>/ { print $1"_"$2"_"$3"_"$4 } /^[A-Z]/ {print $1}' > Alt.fasta

Can anyone help me please?

Answer: Another sed solution:

$ sed '/>/ s/\(.*\)_.*$/\1/g' test.fa

>ART01B_100_M7_ID100005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
Question: Hi,

I am trying to make a Heatmap of some analyzed RNAseq data. I am using ggplot2 with R. The problem I have is that I cannot see well the upregulated and the downregulated genes. I can't get the colors shown to be dark enough. I am doing the heatmap with all the genes of the study organism, around 13,000. I'm also not sure that ggplot2 clusters the genes in my data. Is there a way to cluster the genes and see darker colors in the genes upregulated and the downregulated? The code I am using would be the following:

library(ggplot2)
library(reshape2)

g<-read.csv("datos_prueba3.csv") g3<-melt(g)

plot1 <- ggplot(g3, aes(variable, Gene, fill=value)) + geom_tile() + scale_fill_gradient2(low="dark green", high="dark red", mid="green", midpoint=0) + theme(axis.text.y=element_blank())
plot1
![enter image description here][1]

Best,

[1]: /media/images/f8a7c4c3-4b01-4416-a89d-52e82151

Answer: Doing a heatmap of 13k genes is not going to be informative, it is just too many data. Try the following points:

a) reduce the number of genes to the most meaningful ones, e.g. the differential ones. All others are not different between days and as such do not add any information

b) standardize your data. Expression data have large spans, both if you plot log-scale expression data or fold changes. In this case the extreme values will dominate the color scale (here it is the -20 on the color gradient). Basically the extreme gradient makes it impossible to see differences, e.g. etween -5 and -10. Scaling means Z-scoring. If you have a numeric matrix of log-expression values called `x` then do `t(scale(t(c)))` to get Z-scored data which indicate for every gene the deviation from its mean of all samples.

c) Cluster your heatmap with `hclust` to group similar expression patterns and visualize them. The `ComplexHeatmap` package can do this internally. Alternatively, if you prefer ggplot style use https://github.com/XiaoLuo-boy/ggheatmap which also has parameters to use hclust.

Does that make sense to you?
Question: To make a long story short, I have a bam file that is aligned to the human reference plus a specific contig (human gene + viral vector introduced into a mouse).

The fasta file that was used for alignment got inadvertently deleted, but I have good coverage across the contig. Is there a straightforward way to use the reads aligning to that contig to recreate the contig's fasta file? (All the info should be there!)

Answer: Although reassembly seems the preferred option I thought if Pilon which polishes a reference with standard illumina reads could work . Just for fun I tried following with an ecoli MG1655 bam file.

- Got the length and name of reference chromosome from bam header (in this case only one)
- Created a fake reference with the same name and same number of nucleotides (but all nucleotides were A ; I tried with N and it did not work)
- Then use [pilon][1] and gave the input bam and the fake reference (composed of all A nucleotides) ( I gave as a single end bam but paired should produce better results)

- Pilon generated a fasta (pilon.fasta default name)
- Ran blast2sequence (NCBI;Online) with actual ecoli MG1655 genome)
- Following were the results (99% coverage and 99% identity)
- Dot Plot (https://ibb.co/ifFzjS)

- Alignment details (https://ibb.co/hNCejS)

[1]: https://www.broadinstitute.org/gaag/pilon
[2]: https://ibb.co/ifFzjS
Question: I have a set of cDNA-sequences representing mRNAs and I want to create a gtf/gff file of these. I have tried to use BLAST and BLAT, but I find it difficult to obtain a single reliabe hit that represents all the different exons. And especially with BLAT I get many small hits for each sequence. I also tried to map the sequences using TopHat, but for some reason it crashed, perhaps due to the long sequences.

Any suggestions on how to map such sequences to a genome with reliable information about the exon/intron junctions?

Answer: Use gmap. It was build for that and outputs a gff file with each exon of every sequence in input. Very usefull !
Question: Friends,

I have a microarray dataset AGI IDs in row and samples in column and txt file of transcription factor (AGI IDs), my supervisor asked me to find and extract the expression of only genes existed in my transcription factor list in another word I should extract the rows those are transcription factor. as already I need your help

Thank you

Answer: Take a look at ?intersect.

> v1 <- c("GATA1", "CTCF", "HOXB1")
> v2 <- c("C3", "ATF", "SOX1", "CTCF", "OCT3")
> intersect(v1, v2)
[1] "CTCF"
>

If you are trying to build a subset of a data frame, you can use %in%. See the examples here.
Question: I installed [pygr][1] and its prerequisites.

Running:

```
>>> from pygr import worldbase
>>> dir(worldbase)
['0root', '0version']
```

Where it was supposed to return

['Bio']

according to the documentation (root and version are empty).

Am I doing something wrong?

By the way, this package is five years old and nothing has been contributed to it since. Which leads me to wonder: is the *worldbase* database still maintained at UCLA?

For the curious, what I want to do next is:

```
hg17 = worldbase.Bio.Seq.Genome.HUMAN.hg17()
chr1 = hg17['chr1']
```

[1]: https://code.google.com/p/pygr/

Answer: My experience is that pygr is not maintained anymore. I actually made https://www.biostars.org/p/93364/ as a pygr-compatible interface. You'll find that it works with more or less the same API as pygr, and if you find something that doesn't work the way you want it please tell me and I'll do my best to add it! To work with hg17, 18, or 19 you can just download a multifasta file and pass it to the `Fasta` class. From there you can do what you're attempting and use the reference like a dictionary.
Question: Does anyone know if there is any available code for an interactive web-based ideogram viewer? It would likely be written in javascript. I want to show an ideogram that dynamically displays some annotations (genes related to a phenotype search, for example).

I have found [this ideogram viewer][1], but it is under copyright and the code has not been made available. There is also [Phenogram][2], which looks nice, but is unfortunately only static images. I want to integrate this into a tool I'm building, so existing web-based genome browsers that have ideograms are not sufficient.

Thanks.

[1]: http://bioinformatics.mdanderson.org/ideogramviewer/Ideogram.html?genelist1=TP53&genelist2=ERBB4&mirlist1=hsa-let-7a-2&mirlist2=hsa-let-7c
[2]: http://visualization.ritchielab.psu.edu/phenograms/examples

Answer: I just came across this, and got it up and running quite easily:

https://www.ebi.ac.uk/Tools/biojs/registry/Biojs.Chromosome.html

Pretty much exactly what I was looking for in my original question.
Question: Hi!
We are planning some GBS sequencing experiment. We have the whole genome sequence of organism of interest. But we'd like to test in silico numerous restriction enzymes to choose the best one for our experiment. I wondering, is there any ready program to do that? Particularly, count number of restriction sites in the genome for one or two different enzymes, make distribution of restriction fragment lengths etc. Thanks!

Answer: Apparently, there are several, see these threads:

[Q: In Silico Genome Digestion ?][1]

[Q: Restriction Sites Fragments For Gbs/Ddrad][2]

[Q: how to find the total number of recognition sequence of a restriction enzyme in total human genome ?][3]

Also, see these papers:

https://www.ncbi.nlm.nih.gov/pubmed/27288885

https://www.ncbi.nlm.nih.gov/pubmed/27071032

[1]: https://www.biostars.org/p/948/
[2]: https://www.biostars.org/p/97219/
[3]: https://www.biostars.org/p/173054/
Question: I want to convert a list of fasta ( protein sequences) in a .text file into corresponding nucleotide sequences. A Google search gives me result of DNA to protein conversion but not vice versa. Also, I came across How do I find the nucleotide sequence of a protein using Biopython?, but this is what I am not looking for. Is there any possible way to do it using python.Moreover, I would like to solve it using python programming. I am sure there must be some way to do it rather than writing a code from scratch. Thanks!

Answer: As lieven.sterck points out: this returns you 'a' backtranslation of a peptide sequence. You could use a more dedicated statistical model using codon frequencies from your organism under study, but this is the gist of it:

import random

AA2NA = {
"A": list("GCT,GCC,GCA,GCG".split(",")),
"R": list("CGT,CGC,CGA,CGG,AGA,AGG".split(",")),
"N": list("AAT,AAC".split(",")),
"D": list("GAT,GAC".split(",")),
"C": list("TGT,TGC".split(",")),
"Q": list("CAA,CAG".split(",")),
"E": list("GAA,GAG".split(",")),
"G": list("GGT,GGC,GGA,GGG".split(",")),
"H": list("CAT,CAC".split(",")),
"I": list("ATT,ATC,ATA".split(",")),
"L": list("TTA,TTG,CTT,CTC,CTA,CTG".split(",")),
"K": list("AAA,AAG".split(",")),
"M": list("ATG".split(",")),
"F": list("TTT,TTC".split(",")),
"P": list("CCT,CCC,CCA,CCG".split(",")),
"S": list("TCT,TCC,TCA,TCG,AGT,AGC".split(",")),
"T": list("ACT,ACC,ACA,ACG".split(",")),
"W": list("TGG".split(",")),
"Y": list("TAT,TAC".split(",")),
"V": list("GTT,GTC,GTA,GTG".split(",")),
"*": list("TAA,TGA,TAG".split(","))
}

def aa2na(seq):
na_seq = [random.choice(AA2NA.get(c, ["---"])) for c in seq]
return "".join(na_seq)

print("MARNDCQEGHILKMFPSTWYV*", aa2na("MARNDCQEGHILKMFPSTWYV*"))

One possible output:

MARNDCQEGHILKMFPSTWYV* ATGGCTCGAAATGACTGCCAAGAGGGACACATTCTTAAAATGTTTCCGAGTACCTGGTACGTCTAA

Edit: changed return value of AA2NA.get() for "unknown" amino acids to "---" instead of "-".
Question: Hello,

I would like to have a complete list of genes that are validated targets of microRNA regulation in R.
I have tried getting a complete list of microRNA using the biomaRt package (list of ~1800 microRNA) then input the list to the get_multimir function from the multimir package to retrieve validated miRNA -target gene interaction. In the end I got ~11 000 target genes! This would suggest that over 50% of genes are under microRNA regulation. I guess I did something wrong at some step but I can't figure it out!
Is there a better way to retrieve the complete target gene list?

PS.
This is my first post on this forum, sorry if I didn't pose my question correctly! Any advice is appreciated.

Answer: For the purposes of `multimir` "validated" can include any interaction supported by CLASH/PAR-CLIP/AGO-CLIP data. In this case I'm surprised the number is as low as 50%.
Question: Hi.

I'm using [edgeR][1] for differential expression analysis of RNASeq data (counts downloaded from TCGA).

The workflow I'm using follows the GLM functionality of edgeR, as described in the vignette for edgeR.

```
dge <- DGEList(...)

design <- model.matrix(...) # non-singular

y <- estimateGLMCommonDisp(dge, design)
y <- estimateGLMTrendedDisp(y, design)
y <- estimateTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef = ncol(design) + c(-1,0))
```

Fitting the model, as above, can take ages on our server - seemingly because the dataset is pretty large, the design matrix is pretty large etc etc AND, only a single processor is used by edgeR.

I can't find any options within edgeR that can take advantage of parallel processing and was wondering whether any of the alternative RNASEq expression packages, such as DESeq, had this capability

All the best

Russ

[1]: http://www.bioconductor.org/packages/release/bioc/html/edgeR.html

Answer: I'm currently working on a parallel wrapper for DESeq2. The dispersion estimation for parallel execution is easy, but I need a bit more time to restructure for the LFC calculation.

As Steve mentions, limma is very fast when you have many 100s of samples as in TCGA.
Question: Hi All,

I've tried a few different search engines but I thought I'd put this question out there so there can be an active discussion. Any opinions?

Answer: I would suggest to get a go with MSGF+; it's cross-platform (java) and quite fast.
Question: The GATK output vcf file contains genotypes in which alleles could be separated both by vertical bar ("|") and slash symbols ("/"). For example: 0/1, 0|1, 0/0, 1|1. I found the explanation at [EBI][1] portal. It is said that

> The vertical pipe | indicates that the genotype is phased, and is
> used to indicate which chromosome the alleles are on. If this is a
> slash / rather than a vertical pipe, it means we don’t know which
> chromosome they are on.

Does this statement correct for GATK output?

[1]: https://www.ebi.ac.uk/training/online/course/human-genetic-variation-introduction/exercise-title/want-know-how-we-did-it

Answer: >Does this statement correct for GATK output?

The very page you linked says so. This is part of the [VCF specifications][1], so all correctly formatted vcf files should abide to it:

- /: genotype unphased
- |: genotype phased

[1]: https://samtools.github.io/hts-specs/VCFv4.2.pdf
[2]: https://samtools.github.io/hts-specs/VCFv4.2.pdf
Question: It might be a very stupid question for many of you but, since it's my first variant calling, I didn't figure it out yet.

I have **mpileup**'ped two bam files from two samples, then I filtered the results with **vcfutils.pl** and called the genotypes with **bcftools** **call**. Now I have a **VCF** file containing what I want, but I managed to analyze differences only between the two samples and the reference (which is an assembly coming from a line that is *different* *from both* *samples*).

I have the variants between my two samples and the assembly, but what if I want to detect the differences **between** the two samples? I did it with awk / sed / cut and other command line tools but is there maybe a better and more straight-forward way to do that? Perhaps using bcftools? Up to know, I didn't find it.

Any suggestion appreaciated!

EDIT:
I did it with bcftools gtcheck as well, that works fine. I am asking if there are **other** reliable tools to test!

EDIT 2:
As this post has many views now, I guess it will be useful for everyone to know that I used **bcftools isec** and it worked brilliantly.

Answer: I'm a big fan of vt for this.

http://genome.sph.umich.edu/wiki/Vt#Peek

#summarizes the variants found in sample.vcf

vt peek sample.vcf

Also bedtools is very flexible and can work with VCF files so I would use bedtools intersect for this task.

Others have mentioned database joins to be useful too:
http://www.gettinggeneticsdone.com/2013/11/using-database-joins-to-compare-results.html
Question: Hello, I am trying to cluster a graph but I am running into an error I cannot solve.

I am running the following command:

```
mcl graph.txt -I 1.5 -o outputs/mcl_cluster -te 4 -V all

```

The paths to graph.txt and outputs/mcl_cluster are specified correctly and the graph.txt looks like this:

```
(mclheader
mcltype matrix
dimensions 49347x49347

(mclmatrix
begin

0 4195:0.3275488069414316 $
1 1:0.6507592190889371 $
2 2:0.6507592190889371 $
3 3:0.6507592190889371 $
4 4:0.6507592190889371 $
5 5:0.6507592190889371 $
...
49344 25696:0.3275488069414316 $
49345 47268:0.3275488069414316 $
49346 33405:0.3275488069414316 $
49347 21194:0.3275488069414316 $
)
```
The error message I get is:
```
___ [mclxaRead] (mclmatrix section not found
___ [mcl] no jive
___ [mcl] failed
```
Any idea of why this is happening?

Thank you very much!

Answer: At the very least you are missing a right parenthesis between the `mclheader` and `mclmatrix` sections. It should be like this:

(mclheader
mcltype matrix
dimensions 49347x49347
)
(mclmatrix
begin

That empty line after `begin` may need to be removed as well.
Question: I can't figure out how to use the "-match" syntax even after reading all the documentation I could find. I get these errors:

```
$ cat xml.txt | xtract -pattern Gene-commentary -match Gene-commentary_type:1

Unrecognized argum`ent '-match'
No -element before 'Gene`-commentary_type:1'

$ cat xml.txt | xtract **-element** Gene-commentary -match Gene-commentary_type:1

Unrecognized argument '-match'
No -element before 'Gene-commentary_type:1'
```

**What am I doing wrong?**

---

What I am trying to do is pull the accession of the reference sequence and the coordinates for the region for a given entry in NCBI Gene (see https://www.biostars.org/p/122522/) so that I can run `efetch -format FASTA -seqstart -seqend` and get the appropriate results.

I could parse the XML in python to do it, but it really seems like I should be able to do this in "one line" using entrez direct if only I could get `-match` to work :/

Here is what the XML looks like:

Say I have a gene record in XML

epost -db gene -id 672 | efetch -format xml > xml.txt

According to the outline,

```
cat xml.txt | xtract -outline

1
Reference assembly
RefSeqGene
NG_005905
2

92500
173688
```

I have read:

https://www.biostars.org/p/94539/

http://www.ncbi.nlm.nih.gov/books/NBK179288/ (I followed these instructions to install it, so `which epost` returns `~/edirect/epost`)

http://www.ncbi.nlm.nih.gov/news/02-06-2014-entrez-direct-released/?campaign=facebook-02072014

http://elane.stanford.edu/laneconnex/public/media/documents/EntrezDirect.pdf

Answer: From a bit of experimentation with 'xtract' it appears that the order of the command-line arguments is important, and thus a use of '-match' must be followed by an '-element' option. This appears to be the source of the error message you receive.

Using just 'xtract' the closest I've gotten so far is:

```
cat xml.txt | edirect/xtract \
-pattern Gene-commentary \
-match 'Gene-commentary_type:1' \
-element 'Gene-commentary_accession' Seq-interval
```

You may be able to further anchor the patterns to make the extraction more specific.
Question: Hi all, I am reading paper'A high-resolution map of the three-dimensional' chromatin interactome in human cells. I downlaod their Hi-C data from GEO(GSM1154024) for IMR90 Cell Line.This is used for genome regions interaction.I download the .txt version. The line in this file looks like as follows:

```
HWI-ST216_0375:2:1206:9502:51416#0 chr1 148 - chr19 63789745 -
HWI-ST216_0375:2:1215:19504:87026#0 chr1 6071 - chr13 50507738 -
...
```

As reading above result, I get confused. I know two regions are interacting, like

(chr1 148 -)<->(chr19 63789745-), this two regions should be interacting. But it should be two regions instead of position.I do not know what does 148 and 63789745 means in this dataset.My familiar format should be like (148-500)<->(63789745-63789999),a region map to another region. Any one could help me out how to read result from this paper? Thanks.

Answer: Hi Jesse,

In order to understand this paper, you need to first understand the 3C-based technologies and how the experiment is performed. Please have a look at this review which will help you understand these technologies. You specifically want to focus on 3C part and then move to the HiC the part.

[A decade of 3C technologies: insights into nuclear organization][1]

You also want to read and understand the original HiC paper to understand what a looping interaction means: [Comprehensive mapping of long range interactions reveals folding principles of the human genome][2]

Overview of Hi-C.

(A) Cells are cross-linked with formaldehyde, resulting in covalent links between spatially adjacent chromatin segments (DNA fragments: dark blue, red; Proteins, which can mediate such interactions, are shown in light blue and cyan). Chromatin is digested with a restriction enzyme (here, HindIII; restriction site: dashed line, see inset) and the resulting sticky ends are filled in with nucleotides, one of which is biotinylated (purple dot). Ligation is performed under extremely dilute conditions to create chimeric molecules; the HindIII site is lost and a NheI site is created (inset). DNA is purified and sheared. Biotinylated junctions are isolated with streptavidin beads and identified by paired-end sequencing.

(B) Hi-C produces a genome-wide contact matrix. The submatrix shown here corresponds to intrachromosomal interactions on chromosome 14. Each pixel represents all interactions between a 1Mb locus and another 1Mb locus; intensity corresponds to the total number of reads (0-50). Tick marks appear every 10Mb.

(C, D) We compared the original experiment to a biological repeat using the same restriction enzyme (C, range: 0-50 reads) and to results with a different restriction enzyme (D, range: 0- 100 reads, NcoI).

The topological domains comes later which are basically regions of genome where the elements involved in looping tends to happen in one domain.

I hope this helps a bit.

[1]: http://genesdev.cshlp.org/content/26/1/11.full
[2]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2858594/
Question: I have a file contain multiple gene >sequences. I want to separate into a file with gene ID perl or python script

Answer: As others have mentioned this is answered a lot on the forum. But I can't help myself when it comes to trying to make `bash` do this sort of thing (sequences will have to be linearised).

#!/bin/bash

i=1;
while read line ; do
if [ ${line:0:1} == ">" ] ; then
echo "$line" >> seq"${i}".fasta
else
echo "$line" >> seq"${i}".fasta
((i++))
fi
done < $1

Usage:

$ bash splitfasta.sh multifasta.fasta

Disclaimer:

You should always use a proper parser though (like biopython) as it'll catch many of the special cases. My code just has the bonus of not requiring anything to be installed to run.