cannin/biostars_qa huggingface dataset:
Question: This might be easy but I have not been able to solve it - I have a normal fasta file as below:

```
>seq1
ATCTAC
>seq2
AATCGCATCG
>seq3
ATATACAGC
>seq1
ATCGCGGGC
>seq4
ATTAATTTTAT
```

what I want o do is to reverse the order of lines in the file. So the desired output should be:

```
>seq4
ATTAATTTTAT
>seq1
ATCGCGGGC
>seq3
ATATACAGC
>seq2
AATCGCATCG
>seq1
ATCTAC
```

I tried using ```tail -r file.fa``` but then it does not work because ids get placed after the sequences.

How can I do this?

Thank you

Answer: awk '/^>/ {printf("%s%d\t%s\t",(N>0?"\n":""),N,$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' < input.fasta | sort -t $'\t' -k1,1nr | cut -f 2- | tr "\t" "\n"
Question: Just a very quick question on promoter prediction software. Has anyone used any of the available tools and if so do you have any recommendations on a) the fastest b) the most accurate, and c) the best of both?

Here is a list of some I found! Are any of these scalable to whole-genome data? Are there any others not included here?

Of course if promoter prediction programs aren't the best way to retrieve these data (i.e. promoter sequences +/- first exon), then I'm open to suggestions :) I'm looking to do this for Hymenoptera (i.e. ant, bees, wasps).

Cheers,

Steve

Answer: You may also be interested in looking at some of the following from our labs Links page (http://facility.bioinformatics.ucr.edu/home/links).

Software for Promoter Motif Discovery

Regulatory Sequence Analysis Tools (RSA)

MotifCut (uses Maximum Density Subgraphs)

Verbumculus

DiAlign (aligns sequences with local similarities)

LAGAN Alignment Toolkit Website

FootPrinter

rVista (human & mouse) (dead link)

BioBase (dead link)

ScanACE AlignACE

Consensus, ANN-Spec, Co-Bind, FoldAlign (Gary Stormo Lab)

MEME META-MEME (dead link)

MEME for Arabidopsis (Julian Schroeder Lab)

The Gibbs Motif Sampler Homepage

Gibbs Motif Sampler, Coldspring Harbor

Motif Sampler

Stanford

BioProspector

Collection at Genome Web (dead link)

Link Collection Page from CSL (Michael Zhang)

Promoter Prediction

WWW Signal Scan (dead link)

Neural Network Promoter Prediction

Nucleotide Regulatory Region Analysis

Accuracy (dead link)

Cister: Cis-element Cluster Finder

Human Promoter Prediction: Zhang Lab

Promoter & cis-Regulatory Element Databases

Athena

Eukaryotic Promoter Database

TESS : Transcription Element Search System

Human Promoter Database

Arabidopsis thaliana Promoter Binding Element Database

AthaMap (genome-wide map of putative transcription factor binding sites in Arabidopsis thaliana)

PlantCare (dead link)

PLACE

AGRIS (Arabidosis Gene Regulatory Information Server)

Transfac
Question: Dear all,

after clustering a gene expression dataset (with heatmap.2, or pheatmap), please can you advise, how can i extract the group of gene names corresponding to the dendrogram ? thanks,

-- bogdan

Answer: You can use *cutree* function,

cutree(hc, k)

have a look at following links:

http://stackoverflow.com/questions/20912001/extract-groups-from-distance-matrix-using-cutree-in-r

http://stackoverflow.com/questions/11489696/how-to-know-about-group-information-in-cluster-analysis-hierarchical

hope this will help you !!!
Question: We have a very simple experiment, 1 control (C) and 3 independent treatments (T, B, A). So technically I wrote the following script thinking it would compare C to T, then C to B and finally C to A.

Though at the bottom of my pipeline, it seems that DESeq2 is comparing treatments T vs A. Is there a better way to write the condition so it would compare independentely control to each treatment ? (After that I should cluster genes and output a heatmap)

```
# DESeq1 libraries
library( "DESeq2" )
library("Biobase")
# Heatmap libraries
library( "genefilter" )
library(gplots)

r2 = read.table("matrix.txt", header=TRUE, row.names=1)
head(r2)
samples <- data.frame(row.names=c("C1", "T1", "B1", "A1"), condition=as.factor(c("C1", "T1", "B1", "A1")))
dds <- DESeqDataSetFromMatrix(countData = as.matrix(r2), colData=samples, design=~condition)

# Run DESeq2
dds <- DESeq(dds)

# DESeq2 Results
res <- results(dds)
head (res)

**********

> mcols(res, use.names=TRUE)
DataFrame with 6 rows and 2 columns
type description

baseMean intermediate the base mean over all rows
log2FoldChange results log2 fold change (MAP): condition T1 vs A1
lfcSE results standard error: condition T1 vs A1
stat results Wald statistic: condition T1 vs A1
pvalue results Wald test p-value: condition T1 vs A1
padj results BH adjusted p-values
```

Answer: This has nothing to do with DESeq2. It can't read your mind that you want C1 to be the base level. You just need to `samples$condition <- relevel(samples$condition, "C1")` before making `dds`.
Question: I would like to develop a small scale genomic portal that will (1) accept requests from the user; (2) search several tables of data (currently in the csv format); (3) perform computations in R; (4) output the result in the form of graphs.

This web app will be running on a remote web server (not my local machine) and it should be robust enough that could manage several user requests without breaking... The cBioPortal website http://www.cbioportal.org/public-portal/ is doing something very similar to what I need (but mine is in a much smaller scale)..

I know a lot of Python and R, but I have no idea where to start. I don't know any MySQL but I am willing to give it a go if I need.. Just not sure where to start. Any suggestions?

Answer: I would recommend to consider the following options, which will in the optimal case require very little programming:

- Build your portal based on an existing CMS solution, first choice would be [Drupal][1] + [Tripal][2] depending on the the data in your files, assuming these are genomic data, they could possibly imported easily. The downside is that Drupal is in PHP in case you need to extend it directly. However, you could deploy small standalone python scripts for visualization. And most likely all desired functionality is already covered by extension modules.
- Try to avoid executing heavy computation on the web-server, you need a compute backend of some sort with load balancing (a cluster or cloud solution)
- If R is only for simple plotting I would try to avoid it, and check JavaScript visualization options, there are many good options, e.g. [google chart][3] or [highcharts][4].
- If necessary you could evaluate a R/[Shiny][5] solution for visualization

Note: a beginner cannot implement a secure web-application (especially using DB backend) or set up a hardened server to run on, make sure to get help from a skilled sys-admin and programmer.

[1]: https://www.drupal.org/
[2]: http://tripal.info/
[3]: https://developers.google.com/chart/
[4]: http://www.highcharts.com/
[5]: http://shiny.rstudio.com/
Question: Hi!

I am curious to know, how can I find all the genes within a genomic region, preferable in a file or may be in a separate window.

for instance say, I have a genomic region spanning from chr1:1,783,590-2,597,819 and I want to get all the genes present in this region how can I get it?
I know I can use UCSC browser to see all the genes present in that region but how can I get it as a list in a file.

Thank you

Answer: Or you can create a query in BioMart.
Below an example using the Bioconductor biomaRt package (a wrapper to the on-line web service):

library(biomaRt)
mart <- useMart("ensembl")
mart <- useDataset("hsapiens_gene_ensembl", mart)
attributes <- c("ensembl_gene_id","start_position","end_position","strand","hgnc_symbol","chromosome_name","entrezgene","ucsc","band")
filters <- c("chromosome_name","start","end")
values <- list(chromosome="1",start="1783590",end="2597819")
all.genes <- getBM(attributes=attributes, filters=filters, values=values, mart=mart)
Question: Hello

I have a .ffn file having 1000 sequences. I wanted to check whether all the fasta headers have sequence underneath them, and it would be great if I also get to know about the fasta headers which do not have sequences underneath them.

Thanks

Answer: BioPython is ideal for this.

from Bio import SeqIO

records=list(SeqIO.parse("your_file.fasta","fasta"))
for record in records:
if len(record.seq)==0:
print(record.id)

This will give you ids of headers without any sequences.
Question: I have run the software [sniffles](https://github.com/fritzsedlazeck/Sniffles/wiki) to call structural variants from Pacbio sequences. In the resultant VCF file, most of entries in the ALT field look like:

N[chromosome_3:7199420[

In the example above, this is actually for a variant in a different chromosome, CHR=chromosome_1 and POS=270281, so I guess this is a transponible element coming from chromosome 3 that is present in that location.

I am not familiarized with this format for the ALT field, and was wondering if there is a straightforward way to get the sequence for that element (or any other structural variant found in the VCF). Any ideas?

Answer: What you saw is a BND record, which represents arbitrary rearrangement event with 2 break ends. The `t[p[` format represents “piece extending to the right of p is joined after t” (see details in https://samtools.github.io/hts-specs/VCFv4.2.pdf).

In Sniffles, this is likely a translocation event that might or might not involve transposable element. To extract SV sequence, you can use `-n -1` option to have Sniffles output all SV-supporting reads for each SV entry in the VCF, then you can find read IDs under `RNAMES=`. The first ID usually represents the "primary SV read" that contain representative sequence.

If you are using long read data and are interested in getting all non-reference transposable element sequences based on Sniffles output, you can use TELR (https://github.com/bergmanlab/TELR) that will run Sniffles, find candidate TE loci, and report their sequences based on a local assembly strategy.
Question: I'm interested in extrapolating what transcripts dictates differences between sample in a PCA from RNA-seq data.
[This][1] is the PCA of RNA-seq on different samples (**A**,B,C,**D**,E,**F**,G):

Samples are in duplicate and are circled for better understanding. **Bold** samples are controls, samples connected with arrows are treated samples.

What I want to understand is what are the genes that "moves" **A** to B (Red Arrow) or **A** to C (Blue arrow).
Given the partially similar trajectory of **A** to C and **D** to E the second question would be to check whether some genes in **A** to C are also responsible of the movement of **D** to E.
This same comparison can be made from comparing **A** to B and **F** to G.

Do you happen to know a tool that can allow me to better address this questions?

[1]: https://ibb.co/7Rd81D7

Answer: First of all, there are no trajectories in PCA. It's not a graph/network.

You can get the 'n' list of genes that are most variable between any two groups by rowVars function that is part of several R packages. These highly variable genes are the drivers of PCA. Usually top 500 genes are considered for PCA.
Question: Hello all,
I am working goat genome sample, completed the assembly and masking of genome. Can anyone suggest how to go do gene prediction since augustus doesnt have any ref goat for gene prediction.

Answer: Hi,
1. easy way- use on line training here "http://bioinf.uni-greifswald.de/webaugustus/training/create" and use the paramaeter for you species.

you can prepare training data set from the goat refseq files from NCBI (ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Capra_hircus/latest_assembly_versions/GCF_001704415.1_ARS1/).

2. You can follow the instructions given in the following links for making custom files and training AUGUSTUS.
"https://vcru.wisc.edu/simonlab/bioinformatics/programs/augustus/docs/tutorial2015/training.html"
"http://avrilomics.blogspot.in/2013/04/training-augustus-gene-finding-software.html"
you can use refseq RNA and protein files from the NCBI link for getting RNA and PROTEIN files.

Notes:
I used MAKER to annotate genome initially using RNA and Protein sequences from the closest species and the GFF output from this was used to get 1) multi-exonic 2) non-overlaping 3) relatively non-close sequences (<70% identity) transcripts.
I have converted this GFF file (after filtration) to gb file format using "./scripts/gff2gbSmallDNA.pl" available in AUGUSTUS and used further according to the instructions given in above mentioned links (to be specific first link).
Question: Can somebody help me in converting Gene name to Entrez in R. For example

Gene name: `ANKRD62P1-PARP4P3` has entrez: `23783`

like this page:

http://www.genecards.org/cgi-bin/carddisp.pl?gene=ANKRD62P1-PARP4P3&keywords=ANKRD62P1-PARP4P3

I tried biomaRt but it didn't have the gene (return to NA)!

I tried this:

ensembl=useDataset("hsapiens_gene_ensembl",mart = useMart("ensembl"))

data = "ANKRD62P1-PARP4P3"
ans <- unique(getBM(attributes = c("hgnc_symbol", "entrezgene"),
filters = "hgnc_symbol",
values = data,
mart = ensembl) )

Your help is really appreciated!

Answer: I'd just use org.Hs.eg.db if it's a mapping within EntrezGene that you're interested in.

library(org.Hs.eg.db)
hs <- org.Hs.eg.db
my.symbols <- c("ANKRD62P1-PARP4P3")
select(hs,
keys = my.symbols,
columns = c("ENTREZID", "SYMBOL"),
keytype = "SYMBOL")
# SYMBOL ENTREZID
# 1 ANKRD62P1-PARP4P3 23783
Question: Hi,

I retrieve the Ensembl gene IDs of yeast ribisomal RNA from biomart-Ensembl

```
Ensembl Gene ID
RDN25-1
RDN18-2
RDN5-4
15S_rRNA
RDN37-2
RDN5-6
RDN5-3
RDN58-1
RDN18-1
RDN37-1
```

I need to download the fasta file of these IDs but in Ensembl I could not such a option. Do you have any suggestion please?

Thank you

Answer: You can use BioMart. Follow the sequence below.

Ensembl --> BioMart --> Choose Database --> "Ensembl Genes" --> Select Yeast Genome from list --> Filters (in the left pane) --> Gene --> Input external referenced ID list --> Paste ID's in --> Attributes (left pane) --> Sequence --> Select as needed --> Results Button at top of the page --> Export to "File" as "FASTA".
Question: I want to know whether detected variants are hotspot mutations in TCGA or not.

Though cBioportal provides hotspot site for each gene, it doesn't provide download function.

Are there sites that provide hotspot mutation list or table?

Thank you

Answer: I found two sites providing hotspot information used in cbioportal.
http://cancerhotspots.org/#/home
http://3dhotspots.org/3d/#/home
They are downloadable.
Question: Hello,

I have never made a script in my life.

The ploblem is how to change the fasta names like this input file:

>Glyma04g14800|Glyma04g14800.3
MMLETVAAVPGMVAGMLLHCKSLRRFEHSGGWIKALLEEAENERMHLMTFMEVAKPKWYE
>Glyma05g24460|Glyma05g24460.1
SNVSIDLTKHHVPKNFLDKVAYRTVKLLRIPTDLFFKRRYGCRAMMLETVAAVPGMVGGM

in this output file (change original names to numbers in ascending order, starting with 1):

>1
MMLETVAAVPGMVAGMLLHCKSLRRFEHSGGWIKALLEEAENERMHLMTFMEVAKPKWYE
>2
SNVSIDLTKHHVPKNFLDKVAYRTVKLLRIPTDLFFKRRYGCRAMMLETVAAVPGMVGGM

I'm so grateful for helping.
Regards,
Naka

Answer: awk '/^>/{print ">" ++i; next}{print}' < file.fasta
Question: Hi,
I am a beginner in R. I want to grep "x" from file1 and pipe it with line count command to know how many lines in file1 contain "X". What is the command for it?

Answer: Let's say we have below file:

**file1.txt**

texst1	test1
test2	test1
test3	test1
test4	xtexst1
test5x	test1

We can read the file, and use `grepl` to search for *"x"*. `grepl` returns *TRUE* for matches, then we can `sum` to get the count.

myFile <- readLines("file1.txt")
sum(grepl(pattern = "x", x = myFile))
# [1] 3
Question: I have a list of Uniprot gene IDs associated with Gene Ontology (biological processes), which I have obtained from Uniprot.org. I am showing only one gene ID with associated the biological processes -- because the other genes have a lengthy biological process.

O95831 activation of cysteine-type endopeptidase activity involved in apoptotic process; apoptotic DNA fragmentation; apoptotic process; cell redox homeostasis; chromosome condensation; DNA catabolic process; intrinsic apoptotic signaling pathway in response to endoplasmic reticulum stress; mitochondrial respiratory chain complex I assembly; NAD(P)H oxidase activity; neuron apoptotic process; neuron differentiation; oxidoreductase activity, acting on NAD(P)H; positive regulation of apoptotic process; regulation of apoptotic DNA fragmentation.

**Problem:** Figure out a way to text mining the biological process that is related mitochondria (where mitochondria is mentioned). Would regex be useful to solve this problem? or what other ways that might be useful?

**Expected Result:** the result that I want to get is the following:

O95831 mitochondrial respiratory chain complex I assembly

Your help is appreciated

Answer: +1 for stating your question very clearly.

You might also want to check MitoMiner. It uses UniProt ID/names as the primary ID. It has several templates (or predefined queries). The linked query is "Show all genes and their corresponding proteins for a specific species that are annotated as mitochondrial according to the Gene Ontology (GO)". You need to edit the query to control what information you want to see in the output. Here you can add the GO annotation field.
Question: Hello,
I have a plenty of scaffold and contig names which are non-numeric like Scaffold1 (or any other name. That's why I cannot replace them with number since I have names like Fvb1 etc) instead of "1".

So, can somebody offer a package in R to get Manhattan plots with non-numeric chromosome identifiers? Also, the x axis should be adjustable and not limited to fixed chr number.

Thanks!

Answer: Solved. I used this [GGplot2 Manhattan Plot Function][1]. Thanks to authors!

[1]: https://github.com/pcgoddard/Burchardlab_Tutorials/wiki/GGplot2-Manhattan-Plot-Function
Question: Hello,

A colleague pointed out this problem with GATK, which generates indexes in memory on the fly if it doesn't find those on disk:

INFO 14:05:14,023 RMDTrackBuilder - Creating Tribble index in memory for file dbsnp_132.vcf
WARN 14:09:31,798 FSLockWithShared - WARNING: Unable to lock file dbsnp_132.vcf.idx (could not open read/write file channel)
WARN 14:09:31,798 RMDTrackBuilder - Unable to write to dbsnp_132.vcf.idx for the index file, creating index in memory only

I've read and gone through these two resources:

http://www.broadinstitute.org/gsa/wiki/index.php/Tribble

http://code.google.com/p/tribble/

But I found it a bit awkward having to write my own java class just to generate an index file from a .vcf :-!

VCFTools does not seem to have that functionality either, at least according to their documentation:

http://vcftools.sourceforge.net/docs.html

Anybody knows how to generate those .vcf.idx files in a more straightforward way ?

Thanks in advance !

PD: I just couldn't refrain myself from including this link too :)

http://en.wikipedia.org/wiki/The_Trouble_With_Tribbles

Answer: Depending on what you want to do downstream with these vcf's, you also have the option of using Heng Li's bgzip and tabix. Bgzip compresses and chunks the contents into blocks, then tabix labels the chunks and produces a .tbi index file. Then you can use tabix on the vcf much like you would samtools on an indexed bam.
Question: Greetings all,

Let me start by saying I know how to add read groups to a BAM file. What I am looking for is the easiest way to do it. This problem has been around for a while so I though I would see what neat tricks people are using.

Thanks

Answer: See [this](http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_protocol#Generate_merged_alignment).

The library information of each read group is kept in `rg.txt`

perl -e 'print "@RG\tID:ga\tSM:hs\tLB:ga\tPL:Illumina\n@RG\tID:454\tSM:hs\tLB:454\tPL:454\n"' > rg.txt
samtools merge -rh rg.txt - ga.bam 454.bam | samtools rmdup - - | samtools rmdup -s - aln.bam

Some discussion over [here](http://seqanswers.com/forums/showthread.php?t=4180).

Python: Adding Read Group (@RG) tags to BAM or SAM files:

How can I tell if my BAM file has read group and sample information

Related post: https://www.biostars.org/p/43897/
Question: Hi,

I would like to be able to calculate the relative distribution of introns in genes.

A paper from Jeffares et al, 2006 "The biology of intron gain and loss" shows graph of the proportion of introns and their distribution to show the 5' bias in Eukaryote.

It is explained that : relative position of the intron = bp upstream of the intron in coding sequence / total number of bp in coding sequence. Then, they pool it into ten classes.

I would like to do something similar for at least the human genome. I am interested in doing this for each chromosome to be able to compare the distribution.

I found that the UCSC has a pretty good database with exon position and I was wondering if it was possible to use this data. I know it would probably require to write a program but I am not a programmer and the only code I know is R.

Any help is appreciated,

Thank you

Answer: For a model organism, you could load the TxDb annotation package (see, e.g., TxDb.Hsapiens.UCSC.hg19.knownGene), extract the introns and transcripts from the data, and then calculate the relative offset (of the intron start from the transcript start on the plus strand, for instance)

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
introns <- unlist(intronsByTranscript(
TxDb.Hsapiens.UCSC.hg19.knownGene, use.names=TRUE))
tx <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene,
columns=c("tx_name", "TXCHROM"))

idx <- match(names(introns), tx$tx_name)
offset <- ifelse(strand(introns) == "+",
start(introns) - start(tx)[idx],
end(tx)[idx] - end(introns)) / width(tx)[idx]

The commands are from the GenomicFeatures and GenomicRanges packages. Plot the result

hist(offset)

For less model organisms, a 'TxDb' database can be created with the vignette on the GenomicFeatures landing page.

You could coordinate the offset information with chromosome (with a little bit of trickery to deal with the appropriate factor 'levels' and to focus only on the autosome and sex chromosomes) with

df <- data.frame(offset=offset, chrom=tx$TXCHROM[idx])
lvls <- paste0("chr", c(1:22, "X", "Y"))
df <- df[df$chrom %in% lvls,]

After subsetting factors, it is convenient to drop the unused levels

df$chrom <- factor(df$chrom, lvls)

Visualize the information with, e.g.,

library(lattice)
densityplot(~ offset | chrom, df, plot.points=FALSE)
densityplot(~ offset, group=chrom, df,
plot.points=FALSE)

To aggregate several levels, one possibility is along the lines of

levels(df$chrom) <- c(rep("auto", 22), rep("sex", 2))
densityplot(~offset, group=chrom, df, plot.points=FALSE)
Question: Hi,
Iam currently trying to install EUPAN toolkit (EUkaryotic pangenome analysis toolkit) for pangenome analysis using this installation guide: http://cgm.sjtu.edu.cn/eupan/installation.html
I have downloaded the toolkit but the make command doesn't run and gives error as shown in the screenshot. Kindly guide me step by step commands as to how the issue can be resolved. (This is my first time working in linux)

Answer: I fixed your post.

You are using a version of `gcc` that is too advanced. You have to compile it with a version pre gcc v5.0 (I think) - try gcc v4.1 if possible. I am not sure if you can even install that easily on your Ubuntu 18.04 OS. As you say that this is your first time working in linux, please ask your System Administrator, or a friend / colleague who has experience in this area.

Kevin