cannin/biostars_qa huggingface dataset:
Question: Hi everyone!

I want to perform a gene set enrichment analysis of differently expressed host genes during virus infection. I performed Gene Ontology (GO) annotation with eggnog for the whole host genome and got one .tsv file, an excel file, and two .orthologs files. I've been trying to find literature on how to continue and how to use those files as input for the next tool (can someone recommend a good R package that is compatible with eggnog?), but without success.

I'm very new to this, is there anyone who could give me a hint on how to continue or maybe a protocol/paper that describes the downstream processes after the eggnog annotation?

Please let me know if you need more specific information, any help is highly appreciated!

Thank you!

Answer: I'am new with Bioinformatic so I've been going around in circles for a few months, analyzing my eggnog results
I have a transcriptome and candidate genes I want to do an enrichment test between them, draw graphs, observe the pathway maps and understand the results to which each COG, GO, KO, EC number.... etc. corresponds;

example:

KEGG_ko

kb:K12979
kb:K10603
-
-
kb:K12385

I found an article that talks about TBtools , there is an Eggnog-mapper helper part which is supposed to give separate files for each column with the most important code but unfortunately it does not work with me, maybe that will help you ..
![enter image description here][1]
![enter image description here][2]

[1]: /media/images/aaf66f89-a788-41f7-ba6d-baa9a454
[2]: /media/images/192f0fd1-d578-44a8-b31d-60f7505d
Question: Hello!

Is the source code for the NCBI-BLAST service available anywhere? In particular, I'm interested in producing the visualization NCBI-BLAST offers, such as the one below, using BLAST results from my own installation. I have searched for the source, but can't find it -- it doesn't seem to be included in the BLAST+ tarball, nor can I find it elsewhere on the NCBI FTP site.

Thank you!

Answer: I am reasonably sure that code for the NCBI web interface to BLAST is no longer distributed with BLAST+.

It is still available for the older version of BLAST; the "latest" release (2011) is here: ftp://ftp.ncbi.nih.gov/blast/executables/release/LATEST. Look for the files that begin "wwwblast", but note that the code is not maintained and may be tricky to install.

There are a few other web interfaces to BLAST around (try Google search) which might provide what you want.
Question: The lastz command will look like:

lastz reference.fasta sample.fasta --ambiguous=IUPAC --format=GENERAL > output.fasta

I need to run lastz once for every sample file, because each sample file has to be compared to its own individual reference (I have one reference for each sample).

I'm using find to get all of my 'sample' files, and will do the same for the reference files.

find . -type f -name 'sample.fasta' -not -path "*/temp/*"

The question is, how do I associate each sample file with each reference file, taking as input the output of the find command? Or is there a better way?

Answer: It depends on the location of your reference files. Assuming reference and sample are located within one directory, you won't need associative arrays.

Example file structure:

├── a
│ ├── reference.fasta
│ └── sample.fasta
├── b
│ ├── reference.fasta
│ └── sample.fasta
└── c
├── reference.fasta
└── sample.fasta

Collect paths to samples and references in separate files

find . -type "f" -name "sample.fasta" | sort > samples.list
find . -type "f" -name "reference.fasta" | sort > references.list

Create association sample - reference

paste -d ',' references.list samples.list > file_association.csv

Create commands

while IFS=","
read -r col1 col2
do
x=$(dirname $col2)
echo "lastz $col1 $col2 --ambiguous=IUPAC --format=GENERAL > $x/output.fasta"
done < file_association.csv

or using [GNU parallel](https://www.gnu.org/software/parallel)

cat file_association.csv | parallel --dry-run --colsep ',' "lastz {1} {2} --ambiguous=IUPAC --format=GENERAL > {1//}/output.fasta"

Resulting commands

lastz ./a/reference.fasta ./a/sample.fasta --ambiguous=IUPAC --format=GENERAL > ./a/output.fasta
lastz ./b/reference.fasta ./b/sample.fasta --ambiguous=IUPAC --format=GENERAL > ./b/output.fasta
lastz ./c/reference.fasta ./c/sample.fasta --ambiguous=IUPAC --format=GENERAL > ./c/output.fasta
Question: How can separate one column in two column in R?

id rs143 rs148 rs149 rs1490
1 02003s NA 11 22 11
2 02003s NA 10 11 22
3 02003s NA 11 11 12
4 02003s NA 10 11 11
5 02003s NA 10 11 11

in result i want this format:

id rs143 rs143.1 rs148 rs148.1 rs149 rs149.1 rs1490 rs1490.1
1 02003s NA NA 1 1 2 2 1 1
2 02003s NA NA 1 0 1 1 2 2
3 02003s NA NA 1 1 1 1 1 2
4 02003s NA NA 1 0 1 1 1 1
5 02003s NA NA 1 0 1 1 1 1

Answer: A bit simpler version of Kevin's answer, loop through columns, split, assign column names, and finally bind them back together.

res <- do.call(cbind,
list(df[, 1, drop = FALSE],
lapply(colnames(df[, -1]), function(i){
as.data.frame(list(substring(df[, i], 1, 1),
substring(df[, i], 2, 2)),
col.names = c(i, i))
})
)
)

res

# id rs143 rs143.1 rs148 rs148.1 rs149 rs149.1 rs1490 rs1490.1
# 1 02003s 1 1 2 2 1 1
# 2 02003s 1 0 1 1 2 2
# 3 02003s 1 1 1 1 1 2
# 4 02003s 1 0 1 1 1 1
# 5 02003s 1 0 1 1 1 1
Question: Hello,
I found this code in Biostar handbook (RNA-Seq by Example) in section about genome indexing.

# The reference genome.
IDX=refs/genome.fa

# Build the genome index.
hisat2-build $IDX $IDX

# Index the reference genome with samtools.
samtools faidx $IDX

I was wondering about samtools indexing step. I've never seen anyone using it before. Isn't hisat2-build not enough to properly index genome?

Answer: Yes, the HISAT2 indexing step is enough :)
Question: I am mapping illumina paired-end reads against a transcriptome with HISAT2 and I'm having a moment of confusion understanding one thing which might be easy actually:

Is it normal to have sam flags that point at the same strand for both reads? The technology output should produce one read that maps on the *forward* and one read that maps on the *reverse* strand (when mapping against a genome) but I'm uncertain of what I should observe when mapping against a transcriptome.

I guess the program should map one of the two mates and reverse-complement the other, but does the output flag (and mapping result in general) reflect the mapping strand of the reverse-complement or of the original read? is it normal to have many pairs with both mates on the same strand?

Please help me out of this "theoretical" quicksand.

Answer: It doesn't matter what you're aligning against, you expect PE mates to align with opposite orientations. The strandedness of the underlying data plays absolutely no role in this (your reads don't have strands they have orientations). Unless a very atypical library prep. was used, the flags you should commonly see in the resulting BAM file are 99, 147, 83, and 163. If those don't constitute the overwhelming majority of the flags then something likely went amiss.
Question: I need to run semantic similarity searches on human phenotype ontology.

I am aware of packages like GoSim, SemSim etc. These packages works well with GO. I am looking for a package that can take any .obo file and run semantic similarities on top of it. Do you know of any packages that can do it?

Thanks in advance!

Answer: This is what dnet package can do. See http://supfam.org/dnet/dDAGtermSim.html. However, in this example you need to replace with human phenotype ontology. More information about how to load different ontologies can be found in http://supfam.org/dnet/docs.html.

If you provide details about your annotation information, then probably I can provide more details.
Question: Dear All

I wish to know how can i obtein the Confidence Intervals for Odds Ratio interaction (OR_INT) in epistasis analysis with PLINK.

I'd try with --ci 0.95 but is not work.

plink --file mydata --epistasis --ci 0.95

Thanks in advance!

Adolph

Answer: For future reference, this is answered by https://groups.google.com/d/msg/plink2-users/KMB9YoZqVMQ/7AlJRmobAAAJ .

No, this isn't directly supported by PLINK 1.x (though it's a planned optional column set for PLINK 2.0 --epistasis).

For now, you can infer this from the relation between the beta/odds-ratio and the chi-square statistic:
BETA_INT = ln(OR_INT)
chi-square STAT = (BETA_INT / SE)^2
-> SE = BETA_INT / sqrt(STAT)
95% confidence interval for BETA_INT is [BETA_INT - 1.96 * SE, BETA_INT + 1.96 * SE]; exponentiate these values to get the OR_INT confidence interval
Question: Hi
I'm working on exome dataset I generated variant calls from different combination of aligners vs variant callers. Can i benchmark my exome dataset with Genome-in-bottle (HG001, HG002,HG003, HG004 & HG005)
to find the Truth positives (TP) and False positives (FP). Is this the way ?
Benchmarking mean you can take any query dataset and do comparative analysis with confidence calls (HG001 ....HG005).
Can any one explain, I'm very confused with benchmarking ?

Thanks for your comments.

Answer: You can now benchmark your own variant calls against HG002/HG003/HG004 on [Truwl][1] using 4 popular regions (more coming) and see how your results stack up against the top PrecisionFDA Truth Challenge competitors.

https://medium.com/truwl/accessible-and-uniform-benchmarking-22f598616ef5

[1]: https://truwl.com/workflows/library/Variant%20Benchmarking/0.2
Question: Hello all!
I need to whip up a little network, preferably from known interactions, from a list of about 100 genes in Arabidopsis. All I have are the gene names. I don't care so much about up and down regulation right, just want to represent the genes interacting in one way or another--its more for demonstrative purposes than anything.

Anyone know a quick and dirty way to do this? I don't have any experience in plants so I don't know where to start here.

Answer: You could just go to the [STRING database][1], select "Multiple proteins", paste in your gene names, and select *Arabidopsis thaliana* as organism. That should get you a network in no time.

[1]: http://string-db.org/
Question: I have a sorted, indexed bam of soft-clipped reads. I'm using samtools mpileup and the pileup file it's generating is malformed. For example I see a lot of regions like this, where the # of aligned reads at a position decrements without a corresponding read termination character ($):

```
1 3008463 T 7 .,,.... AGB8A8>
1 3008464 G 6 .,,.$..$ 3GB8<8
```

Oops! I also see plenty of reads that do not have start tokens. Now, granted I'm stress testing it and feeding it really nasty, messy soft-clipped data, so I'm not surprised about the garbage in/garbage out behavior, but if mpileup does not "behave" with soft clipped data I'd like to document this more prominently somewhere, perhaps in the docs for mpileup.

Alright, anybody with any thoughts on this?

-isaac

Answer: The start and end symbols will not appear if the first (last) shown base of a read is filtered out. Have you tried with `mpileup -Q0` to show all the bases? By default, bases with base quality < 13 are not displayed.
Question: I have an `accepted_hits.bam` file which is from tophat output, and I need to use Rsamtools "R package" to do some analysis for the `.bam` file. So any one know how to read this file in R.

Answer: Did you try the package vignette?

http://www.bioconductor.org/packages/2.13/bioc/vignettes/Rsamtools/inst/doc/Rsamtools-Overview.pdf
Question: Hi!

I have calculated the log fold change values for RNA-Seq Data and would like to estimate the significance of the results. I know DESeq does it already, but I want to do it manually after having normalised the counts with RPKM.

Some ideas?

Thanks in advance

Answer: Assuming you've truly done all of the required normalization, then you could just use a T-test or ANOVA (or other applicable linear model). Remember that you'll have lower power than a method like DESeq2 or edgeR since you'll not be using information sharing, but that's the simple manual route.

BTW, why do you want to do this? The various count-based packages are pretty nice and it's usually not a good idea to reinvent the wheel unless you have a good reason.
Question: Im building pipe with secondary structure prediction (spider3). Program do not allow X in fasta sequence. What aa should I use in 2nd structure prediction instead of X? I think about G.

cheers

Answer: What would be the rationale for using G? It has a clear propensity for random coil (more so than strand+helix combined), so there is a chance it will introduce something that is not there. My recommendation would be to use a residue that has almost equal propensity for all three states, and that will hopefully not distort the actual structure in that region. In order, those residues that can be found almost in any three states are: H, R and K. Below are relative propensities for all residues I mentioned.

A better suggestion may be to avoid using sequences with X residues, assuming you have that option. If not, a BLASTp search will tell you what best hits have at those X position, and then you replace them with actual residues from the best match.

coil	strand	helix
G	1.701741723	0.675800457	0.49126143
H	1.12233694	1.0171953	0.896780032
K	0.980548113	0.807582702	1.170839916
R	0.892743354	0.926224811	1.192394115
Question: I know for the reverse procedure, there are lots of methods available. But how can I change my fasta file with sequences like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRQQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRVKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

to like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
QQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
VKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

I know its just about inserting line breaks in each sequence, but is there any way that I can check whether the fasta file is single line fasta or not, and if is, add line breaks in each sequence..

Ps. I am using windows..

Thanks for your consideration

Answer: Here's a quick one in python:

import sys

inputfile = sys.argv[1]
length = int(sys.argv[2])

print '\nConverting ' + inputfile + ' to ' + str(length) + ' chars per line...'

with open(inputfile, 'r') as f:

h = [] #headers list
s = [] #sequence list
for line in f:
if line.startswith(">"): #grab headers
h.append(line.strip().split(">")[1])
else:
s.append(line.strip()) #grab sequences

myseqs = dict(zip(h,s)) #map to dictionary, header:sequence

outfile = open(inputfile.split(".fasta")[0] + '_multi-line.fasta', 'w') #open outfile for writing

def print_seqs(header, sequence):
print >> outfile, '>' + header
while len(sequence) > 0:
print >> outfile, sequence[:length]
sequence = sequence[length:]

for i in myseqs:
print_seqs(i, myseqs[i])
print 'Converted ' + str(len(myseqs)) + ' sequences.'
print 'Done.\n'

outfile.close()

Run as '*python your_input_fasta Nchars'*

**EDIT**

A cleaner, more memory efficient solution, per WouterDeCoster comment:

import sys

inputfile = sys.argv[1]
length = int(sys.argv[2])

outfile = open(inputfile.split(".fasta")[0] + '_multi-line.fasta', 'w') #open outfile for writing

with open(inputfile, 'r') as f:
for line in f:
if line.startswith(">"):
print >> outfile, line.strip()
else:
sequence = line.strip()
while len(sequence) > 0:
print >>outfile, sequence[:length]
sequence = sequence[length:]
Question: I need to do permutation to 10,000 DNA sequences, the length of which are between 200bp~1000bp. In fact, for each sequence, I only need 5,000 permutations to do Monte-Carlo simulation. However, when perl modules List::Permutor or Algorithm::Permute was applied, I found 5,000 permutations of one sequence were quite similar, which went against my intention of random sampling. And to generate a universal set of permutation then randomly sampling is too time consuming. So my question is: is there a way to do a random permutation?

In fact, all I need to do is to search the enriched motif (5mer-10mer) on the sequence sets. So if there is any program for this, plz tell me :)

Thank you!

Answer: How about the EMBOSS program, shuffleseq?

Description: This takes a sequence as input and outputs one or more sequences whose order has been randomly shuffled. No bases or residues are changed, only their order.

You can run it something like:

shuffleseq -sequence myseqs.fa -outseq myshuffledseqs.fa -shuffle 5000
Question: Hello,

I have a large 155mb vcf file that I am trying to sort lexicographically by chromosome number (in the order 1,10, 11, 12 .. etc ) and then by position. I am trying to use vcf-sort from the vcf-tools perl package.

I have been running the command :

perl vcf-sort myfile.vcf > output.vcf

however, the output file loses all chromosome, position, and ref information for some reason and sorts variants by ID. According to the help page, this is suppose to sort the variants in lexicographic order.

I tried the command:

perl -c vcf-sort myfile.vcf > output.vcf

and the output file is sorted in chromosomally chronological order (1,2,3,4..) and all of the necessary information is retained.

Any idea on how to obtain a lexicographical sort by chromosome number (and then position)?

- Parth

Answer: grep '^#' in.vcf > out.vcf && grep -v '^#' in.vcf | LC_ALL=C sort -t $'\t' -k1,1 -k2,2n >> out.vcf
Question: I need a list of polyA site positions. For instance, I extracted TSS positions from a gtf file by taking the starting position of all transcript features. Would the end position of the transcript be the polyA site?

Answer: Many gtf files contain the coordinates for the 3'UTR, too (check the 3rd column). Its end coordinate might be a suitable proxy.
Question: We'd like to perform trio analysis of mice to detect de novo mutations (via WGS and variant calling). The complication: our mouse facility performs breeding in bulk (~15 each males and females). Mom-pup pairs are obviously straightforward. To identify dads, our strategy is to sequence all of the males and infer paternal parentage by comparing dad-pup SNPs. One method is to manually identify SNPs at heterozygous loci for discrimination, but I assume there are tools that can perform this analysis from bulk VCFs. Any suggestions?

Answer: I recently came across [Peddy][1] that I'll be trying out soon. [VCF tools][2] also has a couple of implementations in `--relatedness` and `--relatedness2`, which are based on Yang et al, and Manichaikul et al respectively. There's also [KING][3] which I think was used in ExAC.

[1]: https://github.com/brentp/peddy
[2]: https://vcftools.github.io/man_latest.html
[3]: http://people.virginia.edu/~wc9c/KING/manual.html
Question: I am trying to plot a heatmap of read density around a feature of interest (TSSs) very common in genomics papers. something like this (B):

However, I am struggling a bit in getting to look "right". A bit of background:
I have mapped ChIP-seq reads for pol2 and calculate the coverage, per nucleotide, using bedtools.

coverageBed -d -abam $bamFile -b $TSSs > $coverage.bed
# output:
chr1 67108226 67110226 uc001dct.3 16 + 1 10
chr1 67108226 67110226 uc001dct.3 16 + 2 10
chr1 67108226 67110226 uc001dct.3 16 + 3 10
chr1 67108226 67110226 uc001dct.3 16 + 4 10
chr1 67108226 67110226 uc001dct.3 16 + 5 8
chr1 67108226 67110226 uc001dct.3 16 + 6 8
chr1 67108226 67110226 uc001dct.3 16 + 7 8
chr1 67108226 67110226 uc001dct.3 16 + 8 8
chr1 67108226 67110226 uc001dct.3 16 + 9 8
chr1 67108226 67110226 uc001dct.3 16 + 10 8

Then in R, the genomic position, in column 7, is converted to relative position to the TSS and read counts normalized to the library size. This is converted to a numeric matrix with each row being a TSS and each column the relative nucleotide position. For the plotting the matrix is ordered number of reads per TSS, and the values logged. This is the outcome:

heatmap(cov.mlog, Rowv=NA, Colv=NA, scale="column", labCol = FALSE, labRow = FALSE, col=brewer.pal(9, "Greens"), margins = c(5, 5))

Although the average coverage plot looks as one would expect for pol2, the heatmap is, well, not great. My question is what am I doing wrong and how to improve it? In this and this paper, the coverage is calculate for a bin of nucleotides and not per nucleotide as I am doing. Would that improve visualization? How to do it in bedtools? Should the sorting of the matrix be done differently?

I am aware that could use ngsplot for this, but I am trying to avoid it because my implementation would fit better with my other analysis.

Thank you!

Answer: Hi,

you may find deeptools useful.

Example plots can be found here: http://f1000.com/posters/browse/summary/1094053

Furthermore, deeptools can create bigwig files that are normalized or ratios. As a bonus a Galaxy wrapper is available in the ToolShed.
Question: **Background**

We observed some RNA-seq coverage in regions outside annotated genes - let's called them intergenic. This expression appears to be more pronounced, or unique, to a particular condition.

**Goal**

Find those genomic regions that have higher coverage than expected by random noise alone, along with a read count value (expression). We are looking to identify those regions with high resolution, but rather a broad overview to:

- test whether or not there is a trend for more intergenic expression in some conditions;

- intersect those expressed intergenic regions with other relevant genomic features.

**Data**

Paired-end, total RNA-seq, Vertebrate species, **not human**.

**Possible strategies**

A, näive:

- divide the genome in windows (size?)
- count reads per window
- remove regions containing genes +/- 5kb
- set background: randomly select X regions (1000) with 100 permutations to find distribution of background. Define cut-off as mean (or median + 2*SD).
- Use cut-off to select intergenic regions with high expression. Merge those within 1kb.

B, fancier following a histone mark-style [approach][1]:

- Use `csaw` to calculate coverage using sliding-window (size?)
- remove bins containing genes +/- 5kb
- median coverage across those bins used to filter "expressed regions" (I could also use a permutation approach here)

**Question(s)**

1. Does any of the above options sound reasonable for what I trying to accomplish?
2. Is there some detail missing?

For the window sizes I was thinking about using the average size of exons, since using the size of transcripts could lead to really large windows. Also, if the expression is "transcript-like", short exons - variable length intron - it could lead to large discrepancies in the average coverage and some regions might be missed.

[1]: https://www.bioconductor.org/help/workflows/chipseqDB/#counting-reads-into-windows

Answer: To answer my own question, I ended using two complementary approaches:

1. was to simply count total intergenic reads per condition and compare those (normalizing to mapped reads). The picture as pretty obvious with some conditions showing a lot more intergenic expression than others.

2. I also used `derdinder` to define expressed regions irrespective of annotation, and them used heuristics to select those of interest. We basically had a handful of regions manually annotated, and from those I could derive rules to select other regions of interest. In the end, the rules were stricter than the manual annotation.

Anyway, `derfinder` was a great tool for this purpose, though it was a little fiddly, and there was lot of trial and error with visual inspection on IGV to get the right mix of sensitivity and specificity of expressed regions.