cannin/biostars_qa huggingface dataset:
Question: I have a list of protein GI and would like to get accession number of the genome (DBSOURCE in genbank file) using Eutilities.

What would be the easiest way to get genome accession numbers for a list of protein GI?

Answer: Here is a perl script to do this. This script is explained in detail [here](http://bioinformatics.cvr.ac.uk/blog/?p=276)

#!/usr/bin/perl

use Bio::DB::EUtilities;

#my @ids = qw(817524604 726965494);
my $infile = $ARGV[0];
my @ids;

open (IN,"$infile")||die "can't open $infile\n";

while()
{
chomp($_);
my @ids=$_;
# print @ids."\n";
my $factory = Bio::DB::EUtilities->new(-eutil => 'elink',
-email => 'mymail@foo.bar',
-db => 'nucleotide',
-dbfrom => 'protein',
-correspondence => 1,
-id => \@ids);

# iterate through the LinkSet objects
while (my $ds = $factory->next_LinkSet)
{
#print " Link name: ",$ds->get_link_name,"\n";
my $protid = join(',',$ds->get_submitted_ids);
print "Protein ID:" . $protid ."\t";
#print "Protein ID: ",join(',',$ds->get_submitted_ids),"\t";
my $nucid = join(',',$ds->get_ids);
print "Nuc ID:" . $nucid ."\t";
my $factory = Bio::DB::EUtilities->new(-eutil => 'efetch',
-db => 'nucleotide',
-id => $nucid,
-email => 'mymail@foo.bar',
-rettype => 'acc');
my @accs = split(m{\n},$factory->get_Response->content);
print "Genome Accession: " .join(',',@accs), "\n";
}
}
Question: I am looking for an antibody database where each antibody is linked with its potential/actual target genes/proteins, i.e. a database similar to DrugBank.ca but only for antibodies. It would be really great if the data is in a parse-able format, up-to-date and comprehensive enough to cover most of the common antibodies -- probably asking too much here

Any suggestions on this?

Answer: Hi Arman,

Some of them are (Free+ Paid):

- Antibody Related Databases and Software:

> Antibody related amino acid sequencing tools, nucleotide sequencing tools, structural modeling tools, and hybridoma/cell culture databases can be found below. Speciality research databases that include monoclonal and polyclonal antibodies are also included.

- Human Antibody Initiative
> The mission of the Human Antibody Initiative (HAI) aims to promote and facilitate the use of antibodies for proteomics research. The initiative consists of two separate activities; (1) the generation of a catalogue of validated antibodies from many different sources and (2) a protein atlas for the expression and localization of human proteins in normal and disease tissue. The two separate activities have as their primary deliverables to generate databases with free public accessibility.

> - The Antibody Resource database is aimed to produce a comprehensive catalogue of validated antibodies towards human proteins. This initiative depends on input from a large number of academic groups and commercial companies.
> - The Protein Atlas initiative is aimed to provide comprehensive and annotated database of high-resolution images showing tissue profiles in normal and cancer tissues. Both databases will be open to the public without restriction (no passwords).

- http://www.antibodies-online.com/

I hope this helps a bit.
Question: Hi,

I have few quick question regarding macs14?

- How one can treat the replicates? Do I have to run separately and then take average from the final peaks or is there any way to include it in one run?

- macs14 generates two wig files under "--single-profile" option; one for control and one for treated. which file to use for plotting final peak profile? Do I have to subtract treated - control manually or does it give the final peak wig files?

Thanks in advance

Answer: If you have the replicate, you can take the consistent peaks called in both replicates. There will be variability in ChIP-Seq data, so sometimes you may want to consider union of peaks instead of intersections.

Regarding peak profile, by control do you mean input ? or control in the sense of WT/KO ? You can plot both the treated and control using deeptools2 so that it will be clear that both separated if that's the intention. For UCSC tracks, I would get the bigwigs with -log10(p-value), to show enrichment over inputs.
Question: hi, I would like to create a list of random mm10 genomic coordinates with bedtools. I saw that the useful command to do this is `bedtools random [OPTIONS] -g `. The problem is that they must have different length, similar to my coordinates:

> summary (new $ width)
Min. 1st Qu. Median Mean 3rd Qu. Max.
188 1645 2649 2904 4196 11484

is there any other tool that can set different lengths when generating coordinates?
Thank you!

Answer: There are a couple of different ways to approach this, I think, depending on what you have on hand. If you already have a set of coordinates such as the ones you've entered above, you could simply shuffle these coordinates along the genome using bedtools [shuffle][1]. If you only have a small list of such regions and need a relatively large amount, you could perform the operation repeatedly to get a large number of such regions.

Alternatively it looks like there is an R package for this task. I've not used it myself, but the documentation is [here][2].

[1]: https://bedtools.readthedocs.io/en/latest/content/tools/shuffle.html
[2]: https://rdrr.io/bioc/regioneR/man/createRandomRegions.html
Question: So I have a 6 lists of differentially expressed genes taken at different periods of time. I need a way to show overlap in the genes that recur between lists. Naturally I tried to make a venn diagram, but 5-6 list venn diagrams are hard to follow. Is anyone aware of other methods of presenting data like this? I've been going through google scholars to find papers that present data like this, but I can't find anything.

Thanks, any help is appreciated!

Answer: UpSetR is a recently released R package with powerful graphics for visualizing intersections in complicated cases such as yours.

https://cran.r-project.org/web/packages/UpSetR/index.html
Question: Hi all,

I have an RNA-seq dataset (comparing two different conditions) that I have analyzed using four different methods of differential expression (DE) analysis: DESeq2, limma-voom, liimma-trend, and edgeR. From these 4 different DE analysis, I get the set of raw p values (uncorrected for multiple hypothesis test) for all the genes. I would like to use Fisher's method for p value aggregation to combine the raw p values from these 4 different DE lists and then do the multiple hypothesis correction. Would that be OK? I am wondering because the p values are coming from different programs with different models. I would really appreciate any thoughts or comments. Thanks for the help!

Answer: My two cents / thoughts:

A DE program is supposed to give you a list of DE genes, testing the null hypothesis that gene is not DE. You have your list of DE genes that way. Just pick a DE program and stick with it.

There are a few problems with your approach:

1) Such an approach has not been benchmarked (most DE programs run under default configurations have been benchmarked and shown to work well -- hence they they were published). If combining results from multiple GLM regressions was good and ensured good sensitivity without losing control of the false discovery rate, those programs would already have done it.

2) The Fisher method has [implicitly] the alternative hypothesis that at least one of the programs would call a gene DE (given the info produced by all 4 programs). Is that really what you want to test? You get smaller p-values for each gene with a modest p-value -- more false positives because you're basically increasing the degrees of freedom. I'd use this alternative hypothesis for certain purposes (e.g. aggregating microarray probes for the same gene, or doing a meta-analysis across 10 different low-sample underpowered studies to, in a sense, increase sample size) but not here. In essence, these aggregation methods are used to combine incomplete parts to summarize your results. Using the exact same data AND testing the exact same prediction from the null (i.e. gene counts don't change) for the Fisher method does not really constitute "independence".

3) Related to above, it could be a form of p-value hacking. A t-test might not give me a significant p-value even though my data meets its assumptions so, before running that test, I throw in another test [also one where all the assumptions are satisfied] just to ensure my p-values are lower.

I would therefore not recommend doing this. For similar reasons, I wouldn't recommend making your "DE genes list" by looking at genes with adjusted p < 0.05 in at least one of the 4 programs. You could be more stringent and look at genes with adjusted p < 0.05 in ALL 4 programs, since that's being conservative and those genes have pretty good evidence of being bona fide DE but making a "DE genes list" with a bunch of false positives is much less desirable (I wouldn't even call it a "DE genes list" at that point).

tl;dr I'd recommend that you not do this.
Question: I have a pipeline which looks something like below, its written in bash

Trimmomatic --> Samtools --> Smalt --> Lastz --> custom perl programs.

If anyone of those processes fails they usually return a non-0 exit code, usually a `1`.
However, my pipeline program will continue to run, report the error and return an exit code of 0, indicating success.

How can I get the pipeline to quit/stop with a non-0 exit code when a part of the pipeline fails?

Answer: You're looking for `set -eo pipefail`. Add that after the first (`#!/bin/bash`) line and your script will fail if any part of your script fails.

See for more information: https://vaneyckt.io/posts/safer_bash_scripts_with_set_euxo_pipefail/
Question: I have plant genomic reads (WGS method) through hiseq 2000 with paired-end reads (read1.fastq, read2.fastq). I need to remove chloroplast and mitchondria reads to assemble only nuclear reads. I need to know which is best and fast way to remove contaminants

1. First doing genome assembly and remove mitochondria and chloroplast reads from genome using Blast.

2. Remove first mitochondria and chloroplast reads from fastq files by mapping reads to respective mitochondria and chloroplast read using bwa and get unmapped reads as nuclear reads as fastq file using samtools and picard and then do genome assembly.

Answer: Option 2 is better, mapping with BWA is fast, and you are reducing the total reads to assembly, therefore your assembly time and complexity is reduced.
Question: A colleague of mine and I have just been discussing the meaning of fold changes and though this question has been asked before non of the answers are actually as straight forward as needed so lets try here to solve this once and for all.

I list here what I understand so far and will update it from your answers.

What is the correct way to understand a fold change value of a gene or protein?

- A foldchange describes the difference of two values (eg. difference of expression in gene/protein A between healthy and diseased case)
- Biostatistical porgrams/packages calculate it via: "Log(FC)" = mean(log2(Group1)) - mean(log2(Group2))
- log2 fold changes are used/plotted in graphs as those are nicer to show because they center around 0, giving reductions a negative value and increments a positive value
- log2 fold change values (eg 1 or 2 or 3) can be converted to fold changes by taking 2^1 or 2^2 or 2^3 = 1 or 4 or 8
- To convert the fold change into change in % or anything that is actually tangible/understandable in "real life terms" ... need answers here! (= actual question I want to ask)

Lets get this solved once and for all, im looking forwards to your posts!

Here are great posts explaining more about fold changes:

[graphical representations][1]

[conversion of log2 fold change to fold change][2]

[1]: https://www.biostars.org/p/312980/
[2]: https://www.biostars.org/p/100460/

Answer: Percent change = 2^(fold change) * 100, where the baseline then correspond to 100%

log2FC = 2 : 2^2 * 100 = 400, where the baseline then correspond to 100%

Log2FC = -2 : 2^-2 * 100 = 25, where the baseline then correspond to 100%

Hope this helps
Question: I have list of genes like

LOC109802912

LOC109788332

LOC109789928

LOC109809764
...

and have a fasta file like

>rna-XR_003803819.1 gene=LOC114915746
TGCTTCTTCTCTGTTGGAAGAGTTTAAAAGCAATAAAACTAAGTGTTTTGAGCTCTCTGAAATTGCTGGT
CATGTTGTTGAGTTCAGTGCGGATCAATATGGGAGCCGATTTATTCAGCAAAAGCTTGAAACAGCTActa
ca

>rna-XR_002239520.2 gene=LOC109794983
TCCTATTCATCATGCAGGATGTCAGAATCTTCATTTCAAATTGGACCAATCCTTCTTGTGAATGAACCCT
TCACGATTGATAATGGTCTAATGACACCTACTTTGAAAATTCGAAGAGATAGAGTTGTGGCTCAATACAG
G

how can i get the fasta sequences for genes?

Answer: Try this:

pip install --user pyfaidx

xargs faidx -d " " YOURFASTA.fasta < GENEID_LIST.txt > extracted_sequences.fasta
Question: Hi all,

I am using R script from Getting Genetics Done ( http://www.gettinggeneticsdone.com/2014/05/r-volcano-plots-to-visualize-rnaseq-microarray.html) to generating Volcano plots but my problem is how to show particular single gene name on volcano plot without any condition (i.e foldchange or pvalue). I just want to highlight my gene name "WNT7A" on volcano plot.

It would be great help if someone help me how to modify following code:

> labs=Gene

I tried `labs <- res$Gene[WNT7A]`

but couldnt work...

This is original code from Getting Genetics Done:

res <- read.table("results.txt", header=TRUE)
head(res)

# Make a basic volcano plot
with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2)))

# Add colored points: red if padj<0.05, orange of log2FC>1, green if both)
with(subset(res, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))
with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="orange"))
with(subset(res, padj<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="green"))

# Label points with the textxy function from the calibrate plot
library(calibrate)
with(subset(res, padj<.05 & abs(log2FoldChange)>1), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=.8))

Answer: Perhaps I missed something, but it looks quite obvious to me:

with(subset(res, Gene == "WNT7A"), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=.8))
Question: I have a very fundamental question that I can't seem to find an answer to.

I have a variety of TF and histone marks for a untreated cell line and a treated cell line. I ran a tool that returned a bed file of 'potential' enhancer regions. The untreated cell line identified 36k possible enhancer sites while the treated cell line identified 34k. If I am interested in seeing whether a histone mark or two are enriched / de-enriched at these enhancer sites, should I merge these two datasets together to generate one large dataset? Or should I take only the common enhancer regions (intersectBed) that appear between both datasets?

Why is it 'correct' to merge them together as opposed to identifying common regions? Or vice versa? Or is there something else I should be doing that is entirely different?

I have not worked with untreated vs treated cell types before so I don't wish to proceed too far before determining this.

Answer: I would do check how many enhancer sites are common b/w two conditions

- If high overlap (>90%), I would conclude, treatment has no significant effect and would proceed by the plotting the enrichment of histone marks/TF on the intersection or union (if >95% overlap) of enhancer regions
- If low overlap or very different enhancer sets, I would say treatment has an effect and would elaborate what kind of enhancers are common and what are the "new" enhancer sites plus which ones got lost. You go forward with GO analysis of the neighbouring genes etc, for these groups. Once you are clear with your groups, go forward with the enrichment calculations.
Question: GenCode has a wonderful breakdown of the number of coding/non-coding genes and transcripts here...

http://www.gencodegenes.org/stats.html

Does a similar breakdown exist for RefSeq anywhere? Thanks!

Answer: I don't believe RefSeq do as good a breakdown of transcript types as GENCODE, however, you may be interested in the following resources:

- Comparison
of GENCODE and RefSeq gene annotation and the impact of reference
geneset on variant effect prediction
- NCBI
Reference Sequences (RefSeq): current status, new features and genome
annotation policy
- A
comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in
the context of RNA-seq read mapping and gene quantification

Generally, I think that you'll find that GENCODE is more comprehensive for non-coding genes; however, for the majority of these, exact function is entirely unknown. Most people filter them out of, for example, RNA-seq experiments, in order to (in part) minimise the stringency of a false discovery rate threshold. On the other hand, RefSeq has the feel of a well-curated resource.

Kevin
Question: Hello,

I am trying to figure out how DESeq2 calculates its log2FC measure - we see a strange pattern in the FCs in our data and I'd like to reproduce them by hand from scratch to make sure this pattern does not reflect an error with my pipeline.

Right now, I'm doing this:

deseqOutput<-DESeq(data_collapsedTechnicalReps)
estSizeFactors <- estimateSizeFactors(deseqOutput_DUP)
RLEnormedData <- data.frame(counts(estSizeFactors, normalized=TRUE))
meanOfRLECounts <- data.frame( rowMeans( RLEnormedData[,1:2]) , rowMeans( RLEnormedData[,3:5]) ) # here, condition 1 = cols 1 and 2, condition 2=cols 3,4,5

colnames(meanOfRLECounts)<-c('Condition1','Condition2')

meanOfRLECounts$log2FC <- log2(meanOfRLECounts$Condition1/meanOfRLECounts$Condition2)

Here, condition 2 is the wild type condition - i.e. the samples that I indicate when I use relevel().

Can anyone spot what I'm doing wrong? What data does DESeq2 use to generate its log2FC estimates?

Thank you!

Answer: Take a look at [this thread][1] and the [DESeq2 manuscript][2]

They are using a moderated (shrinkage) estimate.

[1]: http://seqanswers.com/forums/showthread.php?t=49101
[2]: http://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
Question: I'm building a Snakefile and getting an error I'm struggling to fix when trying to add the `log` parameter in my homer annotation rule.

Here is the Snakefile:

threshold="25 33 5 66 75".split()

rule all:
input:
expand("homer_files/annotation_files/ATAC24to31_fetal_hMG_consensusPeak.{THRESH}.homerann.txt", THRESH=threshold),

rule homer_annotation:
input:
"diffbind_files/ATAC24_31_fetalMG_.{THRESH}peakset.bed"
output:
"homer_files/annotation_files/ATAC24to31_fetal_hMG_consensusPeak.{THRESH}.homerann.txt"
log:
"logs/homer/annotation.ATAC24to31_fetal_hMG.consensusPeak.(THRESH}.log"
shell:
"annotatePeaks.pl {input} hg19 > {output}"

Now when I run this without the `log` parameter I get 5 separate files created as expected, however if I try to include it, I get the following error:

SyntaxError:
Not all output, log and benchmark files of rule homer_annotation contain the same wildcards. This is crucial though, in order to avoid that two or more jobs write to the same file.
File "/c8000xd3/big-c1477909/foetal_hMG_analysis/ATAC_24to31_foetal_hMG_May19/Snakefile", line 26, in

I'm not sure what is going on here, and have never came across this before when running very similar rules. The way I'm reading this, a separate `log` file should be created for each iteration of the rule, but the error would suggest otherwise.

Any ideas on how to remedy this would be greatly appreciated.

I'm using Snakemake 5.2.2.

Answer: It was a typo in the braces!

log:
"logs/homer/annotation.ATAC24to31_fetal_hMG.consensusPeak.(THRESH}.log"
Question: Hello,

I am trying to use a for loop in python to loop through a bam file and take out reads that map to certain chromosomes, for further downstream analysis. I am having problems in the beginning, seemingly easiest part of my code.

samfile = pysam.Samfile("mybam.bam", "rb")

positions=["7:151970856", "7:151970856"]

for p in positions:
start = p.split(":")
chr=start[0]
pos=start[1]
chromreads=[]
for read in samfile:
if int(chr) == int(read.tid):
chromreads.append(read)
print(len(chromreads))

Output I'm getting:

2525527
0

Output I expect:

2525527
2525527

It seems like it's not even looping through the second position. What am I doing wrong?

Thank you.

Answer: Answered my own question [from StackOverflow][1]

[1]: http://stackoverflow.com/questions/10255273/iterating-on-a-file-using-python
Question: Deeptools is an incredible package, but I want more plotting versatility.

Anyone kind enough to share workflow/code for plotting computeMatrix heatmaps and average plots in R (ggplot, pheatmap etc.)?

Other suggestions for plotting from .bw coverage files and .bed feature files?

Thanks a lot!!!

Answer: I'm not sure it'd make sense to tidy a matrix of that size for use with ggplot2. Below is how to read a simple file produced by computeMatrix in R (the example file I'm using below can be found in the deepTools source code under `deeptools/test/test_data/`):

m = read.delim("computeMatrixOperations.mat.gz", skip=2, header=F)
m = as.matrix(m[,-c(1:6)])

The resulting matrix is exactly what's plotted in `plotHeatmap`, though you can also use `pheatmap()` if you really want:

pheatmap(m, cluster_rows=F, cluster_cols=F) # This is a useless heatmap, I just use it for testing

Note that there are likes NAs in the data, which cause a lot of things to have issues during clustering.

One big caveat to all of this is that you don't see delineations by region groups and samples. That information is stored in the header:

> h = scan("computeMatrixOperations.mat.gz", n=1, sep="\n", what=character())
Read 1 item
> h
[1] "@{\"verbose\":true,\"scale\":1,\"skip zeros\":false,\"nan after end\":false,\"sort using\":\"mean\",\"unscaled 5 prime\":0,\"body\":1000,\"sample_labels\":[\"SRR648667.forward\",\"SRR648668.forward\",\"SRR648669.forward\",\"SRR648670.forward\",\"SRR648667.reverse\",\"SRR648668.reverse\",\"SRR648669.reverse\",\"SRR648670.reverse\"],\"downstream\":0,\"unscaled 3 prime\":0,\"group_labels\":[\"genes\"],\"bin size\":10,\"upstream\":0,\"group_boundaries\":[0,196],\"sample_boundaries\":[0,100,200,300,400,500,600,700,800],\"missing data as zero\":false,\"ref point\":null,\"min threshold\":null,\"sort regions\":\"no\",\"proc number\":20,\"bin avg type\":\"mean\",\"max threshold\":null}"

That's a json string, which I don't think base R has a function to parse. The important part is that the bit after `sample_boundaries` denotes the beginning and end of each sample (the labels are in `sample_labels`) and `group_boundaries`/`group_labels` does the same for groups of regions. So you can subset matrices in the normal way in R accordingly.

For computing the equivalent of `plotProfile`'s output, you can simply `apply(m, 2, function(x) mean(x, na.rm=T)` (be sure to handle NA values!).

I usually prefer creating plots in R, but if it'd be a good bit of work to just reproduce a moderately complicated equivalent to `plotHeatmap` in R, so I would encourage you to take the advice from https://www.biostars.org/u/16222/ and tweak the python code in deepTools instead.
Question: I got the good quality genome IDs for 54000 genomes like below:

#genome
G001281285
G000014725
G000775715
G000254175
G001380675
G900057405
G001076295

and I also got all 74000 genome sequence files compressed in fna folder like below

cd fna/

G001284865.fna.bz2 G002910165.fna.bz2 G009390615.fna.bz2
G001284885.fna.bz2 G002910195.fna.bz2 G009390655.fna.bz2

Now could you please help how I can extract the 54000 genome sequence files based on above genome IDs from fna/ folder?

Answer: for i in $(cat list.txt); do mv "$i".fna.bz2 fna/filtered/; done

Where list.txt is your list of high quality genomes and filtered/ is a new directory.