cannin/biostars_qa huggingface dataset:
Question: Hi,

I'm trying to run the PathSeq software and I followed all the steps exactly, but I'm getting the following message:

Client.AuthFailure: Not Authorized for images: ami-5fb67036

I'm new to cloud computing and I have no idea what the problem is. Any help would be greatly appreciated!!!

Thank you!!!!!!

Anna

Answer: The problem was with the image that the Broad Institute provided. They have fixed the problem.
Question: I have two: Table-A and Table-B. I want to extract the contents of Table-B matching Table-A something like Table-C for all the Motif_ID of column-1 of Table-A. Can anybody help me with R-script.

Table-A

```
Motif_ID Size Overlap Overlap Fold enrichment P value FDR Motif_ID DBID TF_Name Family_Name
M5301_1.01 181 181 14 16.19 6.81E-18 6.12E-15 M5301_1.01 ENSGALG00000013342 GBX2_CHICK Homeodomain
M4396_1.01 500 500 16 6.7 1.53E-14 4.57E-12 M4396_1.01 ENSGALG00000005048 SMARCC1 Myb/SANT
M4435_1.01 493 493 16 6.79 1.20E-14 5.38E-12 M4435_1.01 ENSGALG00000010036 FOSL2_CHICK bZIP
```

Table-B

```
Motif Gene Start Stop Log-odds p-value Site
M0082_1.01 Chr25_scale5_25_989144_990145_+_-R 380 389 - 11.5048 2.30E-05 AGCCTCAGGG
M0082_1.01 Chr25_scale6_25_989242_990243_+ 515 524 + 11.5048 2.30E-05 AGCCTCAGGG
M0082_1.01 Chr25_scale9_25_997105_998106_+_-R 378 387 - 11.5048 2.30E-05 AGCCTCAGGG
```

Table-C

```
Gene Start Stop Log-odds p-value Site Size Overlap Overlap Fold enrichment P value FDR Motif_ID DBID TF_Name Family_Name
M5301_1.01 Chr25_scale1_25_981718_982719_+_-R 516 529 + 10.3933 8.65E-05 TAATTTGCTGATTA 181 181 14 16.19 6.81E-18 6.12E-15 M5301_1.01 ENSGALG00000013342 GBX2_CHICK Homeodomain
Chr25_scale2_25_981736_982737_+ 455 468 - 10.3933 8.65E-05 TAATTTGCTGATTA
Chr25_scale3_25_985474_986475_+_-R 758 771 + 12.5796 2.05E-05 TAATTTGCCCATTA
Chr25_scale3_25_985474_986475_+_-R 758 771 - 13.426 1.11E-05 TAATGGGCAAATTA
Chr25_scale4_25_985197_986198_+ 508 521 + 13.426 1.11E-05 TAATGGGCAAATTA
Chr25_scale4_25_985197_986198_+ 508 521 - 12.5796 2.05E-05 TAATTTGCCCATTA
Chr25_scale5_25_989144_990145_+_-R 523 536 + 12.5796 2.05E-05 TAATTTGCCCATTA
Chr25_scale5_25_989144_990145_+_-R 523 536 - 13.426 1.11E-05 TAATGGGCAAATTA
```

Answer: merge(Table-A, Table-B, by.x="Motif_ID", by.y="Motif")

This should work, however IDs in `Table-A$Motif_ID` **must be unique**.

PS., similar questions were asked [here][1] before. Also I don't know how merge will cope with duplicate column names (get rid of one `Motif_ID`).

[1]: https://www.biostars.org/p/91074/
Question: Hello community,

I am using snakemake to make a pipeline. I want to add the bowtie2-build from bowtie2 to my current snakefile as follow:

rule bowtie2Build:
input: "refgenome/infected_consensus.fasta"
output: "output/reference"
shell: "bowtie2-build {input} {output}"

So I should be expecting the following files:
reference.1.bt2
reference.2.bt2
reference.3.bt2
reference.4.bt2
reference.rev.1.bt2
reference.rev.2.bt2

But it seems the problem lies in the output. How can I write the output?

Answer: You wrote this:

rule bowtie2Build:
input: "refgenome/infected_consensus.fasta"
output: "output/reference"
shell: "bowtie2-build {input} {output}"

What snakemake is going to do is to check if the file output ( in this case "output/reference" ) exist after executing the rule. Which doesn't since it is the basename for bowtie

What you can do is to pass the basename for the index as a parameter to snakemake function. Something like the following:

rule bowtie2Build:
input: "refgenome/infected_consensus.fasta"
params:
basename="output/reference"
output:
output1="output/reference.1.bt2",
output2="output/reference.2.bt2",
output3="output/reference.3.bt2",
output4="output/reference.4.bt2",
outputrev1="output/reference.rev1.bt2",
outputrev2="output/reference.rev2.bt2"
shell: "bowtie2-build {input} {params.basename}"
Question: Hello!

I have MBD-seq datasets of Dnmt1 Knock-Out and control cells.

As expected, Dnmt1 KO covered very few genomic regions compared to control sample since no Dnmt1 presented in the KO. The problem is that the reads in the KO sample were concentrated in the "few genomic region" made the intensity too high.

What I'm curious here is how I should normalize such data which the samples are expected to be different overall reads?

For example, if KO and control have 100 and 10,000 detected bins (> 0 reads), and each of them has a million total reads, each bin in KO will 100 times higher reads leading to the biased quantification of MBD-seq enrichment.

Would it be OK if I subset 1/100 of the reads from a KO sample to compensate the differences?

How do everyone think?

Thank you!

Answer: I suggest you go through the manual of MEDIPS https://bioconductor.org/packages/release/bioc/html/MEDIPS.html to get a guide line for the analysis. It also covers normalization.
Question: Hello!!

I have been trying to figure out how can I calculate the number of SNPs using sliding windows?
I have a list with three columns: Scaffold"\t"Scaffold_Length"\t"Number_SNPs (per scaffold)

Scaffold_28	70818817	894731
Scaffold_3	5123947	57985
Scaffold_13	4491039	67622
Scaffold_12	3793473	51663
Scaffold_23	3593776	31841
Scaffold_11	3547442	63973
Scaffold_26	2720936	36018
Scaffold_16	2719413	24318
Scaffold_27	1987753	53938
Scaffold_24	1647859	18408
Scaffold_9	1630703	15792
Scaffold_32	1545880	21094

. . .

Based on the second column, I want to use sliding windows (500kbp) and calculate how many SNPs are there into each sliding windows.
I performed: bedtools makewindows, but I have not figured out how to count and sum the SNP density.

Thanks a lots for your help!!

Answer: One way is to pipe BED-formatted, sliding windows into `bedmap`, using its `--count` operator to count the number of SNPs that fall within each window.

The following `bedops` statement would generate 500knt windows from scaffolds, spaced every 100knt. These windows are passed along to `bedmap`, which counts the number of SNPs that fall in each of those windows:

$ bedops --chop 500000 --stagger 100000 -x <(awk -vOFS="\t" '{ print $1, $2-1, $2; }' scaffolds.txt | sort-bed -) | bedmap --echo --count --delim '\t' - <(vcf2bed < snps.vcf) > answer.bed

The result is written to `answer.bed`, each line of which containing the window and the number of SNPs over that window.
Question: I have stranded RNA-seq data which doesn't look stranded when I visualize bedGraphs in Genome Browser:

![enter image description here][1]

I tried different aligners (star, bowtie, bwa), different ways to make bedGraph files, nothing works. But when I run salmon on these samples, it detects ISR (paired-end stranded library). I'm very confused why it doesn't look stranded in the browser but I know is stranded. Any help will be really appreciated! I've been trying to solve it on my own for a while...

[1]: https://s8.postimg.cc/4yxzbs4bp/strand_error.png

Answer: Hi,

The reason for this is that whether or not to consider a fragment for quantification in salmon is considered *after* the mapping is written. That is, the reads are mapped to the transcriptome, without regard for the inferred or provided library type, and then these mappings are written to file if requested by the user. However, when salmon considers the probability of these mappings arising from different transcripts, it will discard those that are not compatible with the library type. So, my guess is that you're seeing in this output multi-mapping reads where the sequence matches, and hence there is a possible allocation of this read to the transcript, but for which these mappings are later discarded during quantification. You could, for example, filter the mapping file based on the tags of the alignment records, to only keep those reads that align according to the ISR library type.
Question: Hi all,

I have BLAST+ version 2.11.0+ on my Linux computer and I would like to update this version to the latest (2.13.0+).

I downloaded successfully the newest version from the [NCBI website][1].

However, after carefull reading of the manual, I'm still having trouble configuring this new version instead of the previous one..

When I type from anywhere in a terminal:

```
$ blastn -version
$ blastn: 2.11.0+
```

whereas when I go directly into the bin folder of what I've downloaded:

```
$ blastn -version
$ blastn: 2.13.0+
```

**So, the issue is with the PATH variable.** To try to solve it, I did:

$ export PATH=$PATH:$HOME/ncbi-blast-2.13.0+

but it still doesn't work with the new version..

I know I always have trouble with the PATH variable when I'm installing/updating software. So if you have an idea on how to solve that, I would really appreciate it!

Many thanks for your help,

[1]: https://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/

Answer: You should prepend instead of appending the new directory to your PATH:

export PATH=$HOME/ncbi-blast-2.13.0+:$PATH

Now, when you run your blastn command, the shell will look in the 2.13.0 folder first.
Question: Hi all,

I get the following Error:

Error in do_one(nmeth) : NA/NaN/Inf in foreign function call (arg 1)

I'm trying to make a heatmap (pheatmap package) using K means.
I have a large matrix (x) and I did the following:

x[x < 0] <- NA

x2 <- x[complete.cases(x),]

RedBlue <- colorRampPalette(c("blue", "white", "red")) (20)

pheatmap(x2, scale="row", kmeans_k=100, cluster_rows=TRUE, cluster_cols=TRUE, clustering_method = "ward.D2",color=RedBlue, legend=FALSE, show_colnames=TRUE, fontsize=15, fontface="bold", border_color=NA, width=20, heigth=100, cellwidth=30, cellheight = 5, show_rownames = FALSE, filename = "Heatmap.tiff")

I checked for NA/NaN/Inf but I don't have this is my matrix. What can be another reason for this error?

Thanks!

Answer: This is not a bioinformatics question. Also before posting in a more appropriate channel, provide a reproducible example.
Now since I am here...
This error can happen for several reasons the most common being:
- presence of NAs, including those produced by scale() in the case of variables with 0 variance
- presence of non-numerical values
Question: I want to download the annotation file in gff3 format for the corresponding genome. As this fairly easy on the ncbi-webpage I don't find a possibility to do the same with efetch or the kind.

I hoped I could use something like this:

esearch -db nuccore -query "$genome_id" | efetch -format gff3 > "$path_data/$genome_id.gff"

Answer: There are a couple of strategies you can try, depending on what you mean by $genome_id. In each case, it's a matter of finding the right FTP path, and then using wget to get the *genomic.gff.gz file in that path:

1. If you have assembly accessions, you can get FTP paths for each from
the assembly_summary.txt file, and loop through them with wget. See
https://www.biostars.org/p/61081/ for a good post on the approach
2. If you have nucleotide sequence accessions for chromosomes, you can
use esearch to directly query the Assembly database, and get the FTP
path from the document summary:

esearch -db assembly -query NC_000913.3 | esummary | xtract -pattern DocumentSummary -element Taxid,Organism,AssemblyAccession,FtpPath_RefSeq

3. If you have nucleotide sequence accessions that don't directly work for queries in the Assembly database (e.g. contigs or scaffolds), you can query in nucleotide first and link to assembly:

esearch -db nuccore -query NZ_GL379776.1 | elink -target assembly | esummary | xtract -pattern DocumentSummary -element Taxid,Organism,AssemblyAccession,FtpPath_RefSeq
Question: I am working on NCI-60 database for data mining purpose. In NCBI I found miRNA expression profiles ([link][1]). But for mRNA of each cell lines I found 2 kinds of expression profile array in NCBI ([link][2]) categorized as Affymetrix HG-U133A and HG-U133B. They also overlap many mRNA probes and their expression values. I can't understand which one to use or how? What is the difference?

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26375
[2]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE5720

Answer: [This PDF][1] from Affymetrix clarifies the difference between 133A, 133 B and 133 plus. As for using the datasets, look at this question: https://www.biostars.org/p/56657/. It seems 60 samples are analyzed using both the chips.

Copy/pasted from Affymetrix website:

> The HG-U133A Array includes representation of the RefSeq database sequences and probe sets related to sequences previously represented on the Human Genome U95Av2 Array. The HG-U133B Array contains primarily probe sets representing EST clusters.

[1]: http://media.affymetrix.com/support/technical/datasheets/human_datasheet.pdf
Question: I've assembled a de novo transcriptome from RNAseq data. I'm comparing it against a reference cDNA set that is publicly available. However, I'd like to see if my assembly has produced any novel transcripts. Is there something like a reverse alignment? I'd like to find transcripts that I've assembled that are not included in this file of cDNA. I will then map them to the genome using blat or a similar tool.

How can I do what seems like the inverse of an alignment?

Answer: Check here :

https://www.biostars.org/p/14863/
Question: Hi,

I have a set of transcripts that a collegaue has identified as being up regulated in a study species cf. a reference species. Each transcript has a set of GO terms associated with it. The GO Terms were pulled from a gb file.

The data is in R in a list, names are transcript ID, and the values are the GO terms.

e.g.

$WS010775.1
[1] "0042567" "0005634" "0007155"

$WS007996.2
[1] "0005834" "0005886" "0003924" "0004871" "0010659" "0071456" "0042462" "0007186" "0007602" "0008104" "0007165" "0007268"

$WS010604.4
[1] "0016021" "0007155"

I'd like to present a summary of biological function/features that are represented by those up regulated transcripts. I found REVIGO (http://revigo.irb.hr/) but I believe that it requires that I do GO Term enrichment first to get the p-values that REVIGO requires. Are there any tools I can use to do this, given that the species I'm study is not human/mouse/etc? Tools I've looked at find the terms and do enrichment analysis. But I already know the GO terms of each transcript, I just want to calculate enrichment.

Thanks,

Ben W.

Answer: One of the easiest GO enrichment tests to apply is the Fisher Exact Test approach. Take a look at this response: http://stats.stackexchange.com/a/72556/5143. This is easily done in R.
Question: This question might not be exactly on-topic here, but I thought probably most of you have been in a similar situation before. When describing your methods in a manuscript, do you cite the components of a larger script? For example, MITObim uses MIRA to do mapping, so should I cite MIRA if I used MITObim in my work, or would citing MITObim be sufficient? MITObim is more than a simple wrapper script, but it relies heavily on MIRA.

Just wanted to know what etiquette says in this case; I haven't written a real paper before.

Answer: If it is a published tool, cite only that tool and not anything that tool is. The publication of the tool will used should have citations for anything it used. In some cases you can adjust settings to programs inside that tool, you should make sure to mention any non-default settings.

However, some tools allow you to choose which sub-program you use. For example RSEM can use STAR, bowtie or bowtie2. You should cite this because it allows readers to see exactly what software you used to perform your analysis. RSEM can also take SAM/BAM files as direct input, so if you used something totally different, you should obviously cite it. In your case, if you selected MIRA over another program for MITObim, you should cite it. If MITObim can use MIRA and only MIRA, you probably shouldn't cite it.

If it is a custom script, cite the tools or libraries involved. As well as the author

A fuzzier case may be in instances where the item in question is more protocol than software (e.g. Trinotate). In that case you should probably cite the tools used.

There's always a judgement call in what level you stop at. The important thing is to generally cite things that contributed specifically to your analysis. I guess one way to think of it is what programs did you (or other authors) personally use directly, either by running the program or by specifying it specifically over another program. Giving credit is obviously important, but when dealing with methodological items (wet or dry) the most important thing is giving the reader enough information they could repeat the analysis.
Question: How can I stack several manhattan plots in the same plot looks like this https://imgur.com/a/pj40c using R?

That plot is made using excel and not follow their real map position.

In my case I have 4 p-values from 4 different phenotypes. They are estimated from same SNPs. I want to stack them in the same plot but distinguish phenotypes by different colors.

# Input data format, example.
SNP	CHR	POS	pval_1	pval_2	pval_3	pval_4
a	1	100	0.1	0.5	0.2	0.1
b	2	110	0.2	0.6	0.3	0.5
c	3	120	0.3	0.7	0.1	0.6
d	4	130	0.4	0.1	0.4	0.2

qqman is a very useful R package, unfortunately it is only allowed for single trait.

Answer: With a bit of customization you could do it in R with ggplot:

library(dplyr)
library(ggplot2)

tmp <- data_frame(
POS=seq(1,10000),
pval1=1/runif(n = 10000, min=0, max=0.2),
pval2=1/c(runif(n = 9999, min=0, max=0.2), 1e-5),
pval3=1/runif(n = 10000, min=0, max=0.3),
pval4=1/runif(n = 10000, min=0, max=0.4))

tmp.tidy <- tmp %>% gather(key, value, -POS)

ggplot(tmp.tidy, aes(POS, value, color=key)) +
geom_point()

**Example**:

[https://ibb.co/mYnb3w][1]

Use the + facet_wrap(~CHR) to create a graph per chromosome.

[1]: https://ibb.co/mYnb3w
Question: Common Workflow Language (**CWL**) https://github.com/common-workflow-language/common-workflow-language / http://common-workflow-language.github.io/draft-3/ has been trending on my twitter timeline during the last weeks.

However the spec is quite large and I find it hard to get some simple examples.

Furthermore, I have the feeling that all engines require a lot of dependencies or docker. I'd like to test my makefile-based workflows using CWL, how should I write and test the following simple **Makefile** using CWL:

```
SHELL=/bin/bash
.PHONY: all
all : database.dna

database.dna : seq1.dna seq2.dna seq3.dna
cat seq1.dna seq2.dna seq3.dna > database.dna

seq3.dna : seq3.rna
tr "U" "T" < seq3.rna > seq3.dna

seq3.rna :
echo "AUGCGAUCGAUCG" > seq3.rna

seq2.dna : seq2.rna
tr "U" "T" < seq2.rna > seq2.dna

seq2.rna :
echo "AUGAAGACUGCGAUCGAUCG" > seq2.rna

seq1.dna : seq1.rna
tr "U" "T" < seq1.rna > seq1.dna

seq1.rna :
echo "AUGAAGACUGACUCGUCG" > seq1.rna
```

**EDIT**: feel free to add the file for your favorite workflow-engine as an answer.

https://twitter.com/PaoloDiTommaso/status/625995681434607616

https://twitter.com/smllmp/status/625999447869231104

Answer: Responded with an example on this github issue: https://github.com/common-workflow-language/workflows/issues/1
Question: Hi All,
I intersected to data sets and I have two columns like:

m01 m01
m01 m02
m01 m05
m01 m032
m01 m02
m01 m01
m02 m06
m02 m01
m02 m02
m02 m09
... ...
m0500 m023
...

I would like to get number of matches with others like:

m01 m02 m03 ......

m01 25 45 98 .....
m02 90 223 12 ......
.
.

Would you please help me how I can do that?

Thank you very much

Answer: You didn't say which language you are using. I guess since you tagged it with awk and perl, that's what you want?

Anyway, if you were using R, you could do it with a call to table:

## Maybe you have this in a tab delimited file already
> dat <- read.table('/path/to/2-column-file.txt')
## but I'll generate a random table that looks close-enough to your data:
> set.seed(1)
> dat <- data.frame(x1=rep(c('m01', 'm02', 'm03'), each=10),
x2=sample(c('m01', 'm02', 'm03'), 30, replace=TRUE))
> head(dat)
x1 x2
1| m01 m01
2| m01 m02
3| m01 m02
4| m01 m03
5| m01 m01
6| m01 m03

> table(dat$x1, dat$x2)
m01 m02 m03
m01 3 4 3
m02 2 3 5
m03 4 4 2
Question: Hi, I have one BAM file which contains all alignments (include those not used in variant calling, such as non-PF, non-mapping and duplicate reads) generated for an assembly. How to filter these useless mapping? I know that Picard MarkDuplicates can be used to remove duplicates.Thank you.

Answer: The bamtools package offers a wide range of filters, including user-definable filters defined in JSON notation. It includes filters for reads failing vendor QC, unmapped reads and pre-marked duplicates.
Question: This is a general question about results that I have seen several times but only recently considered. If you do something like ChIP-Seq or HITS-CLIP or other techniques that immunoprecipitate DNA or RNA bound by a protein, then run motif analysis on the areas of the genome mapped to the reads that are pulled down, but do not find the binding motif for the protein that you used for the IP, what does it mean? Is this expected? Did something go wrong with the wet lab steps, or analysis?

Answer: In my opinion it suggests one of two things:
1. Either your antibody has specificity issues. It is not uncommon that antibodies are pulling down other areas than the expected. TSS in particular. Best control for that is to look for signal in a knockout cell line. I know some people who raised five (as I remember it) antibodies against a target and got quite consistent peaksets in chipseq. They later discovered that all antibodies resulted in plenty of peaks in a cell line that did not have the gene encoding the target, suggesting that the peaks were due to some crossreactivity or 'sticky' regions. Ouch... Not a good day in lab.
2. Alternatively the antibody recognizes the target properly, but the assumed motif is not right or the majority of binding is due to indirect recruitment. Again, data from a knock out cell line would be handy to convince your audience of that.
Question: I've got the DNA mutation information which can be in a file format like VCF or MAF. Actually what I want to look at is how the mutation influence the corresponding amino acid changes.

Are there any tools can generate the mutated protein sequences corresponding to the mutated DNA in batch?

Many thanks.

Answer: Great question (for that reason I upvoted it). You might try the R package "CustomProDB" which takes a VCF as input and outputs a protein fasta incorporating SNVs and Indels.

http://www.bioconductor.org/packages/release/bioc/html/customProDB.html
Question: I couldn't find much in the SRA hand book about the differences between the read prefixes. I am downloading the Tardigrade reads and was just curious (https://www.ncbi.nlm.nih.gov/sra/DRX012456[accn])

Answer: They were submitted to different databases. SRR to NCBI, ERR to EBI and DRR to DDBJ, I believe.