cannin/biostars_qa huggingface dataset:
Question: Hi All,

I would just like some clarification of terminology regarding a detail of gene coexpression network construction. Let's say I have two RNA-seq datasets, each dataset containing `n` replicates, and each dataset representing sequencing data from the same biological system in two different experimental conditions. How should I construct the data matrix for input to something like WGCNA if I want to analyze gene coexpression networks ***across experimental conditions/interventions?***

What I imagine is that each row of the matrix represents data from one gene, and each column represents data collected from one of the replicates in an experimental condition. So for example, one particular row of the matrix would look like this:

c1R1 ... c1Rn c2R1 ... c2Rn
gene x [val, ... val, val, ... val]

Where the first column `c1R1` corresponds to the data from the first experimental replicate in the first condition, and the last column `c2Rn` corresponds to the nth experimental replicate in the 2nd experimental condition. For coexpression analysis, each row is then correlated with every other row in a pairwise fashion, an adjacency matrix is constructed from the correlation analysis and then other analyses such as module detection can be conducted based on the resulting adjacency matrix.

I just want to verify that this is an appropriate method for organizing data if one wishes to construct coexpression networks for genes "across an intervention".

Answer: Hi mantale1,

That's exactly correct. By including replicates from both conditions, network will reflect both the specific pathways that are co-regulated during your condition of interest, as well as whatever genes are constitutively expressed in the organism.

If you were to then start added samples from other unrelated conditions, you would be both improving the accuracy of the global co-expression network due to the increased information, but also would be reducing the signal resulting from the intervention you are interested in.

Couple things you might consider:

1) Depending on the number of replicates you have for each condition, you may end up with a very noisy co-expression network. Most of the methods were developed for microarray data where you are likely to have many more samples. With less then 10 replicates across both conditions, you are likely have a large number of spurious correlations.

2) You might consider filtering out genes which are not differentially expressed across your intervention. This will help both with eliminating spurious correlations, and also help to bring out the signal specifically due to the intervention.

Keith
Question: Hi,
I have 3 metagenomes, all of which are from enrichment cultures. My aim is to assemble the genomes of bacterial strains I have not been able to isolate which may degrade my compound of interest. I thought it would be best to merge the samples beforehand and this would improve the coverage of the MAGs but after reading some other posts I'm not so sure. Would it be best to concatenate samples? If so, is it ok that I have done this prior to the trimming stage, i.e I have concatenated all the r1.fq.gz and all the r2.fq.gz files and now plan to trim these files?
Thanks,
Jess

Answer: Hi Jess and welcome at Biostars.

To make it short: Yes your approach very likely will work. I do similar merging operations when I get WGS metagenome data from, for example, a time series experiment.

To outweigh the posts you mentioned, maybe take a look at [metabat][1] or [concoct][2], two software packages that can be used to isolate genomes from you metagenome. Though both essentially start binning using k-mer frequencies, they also assume bam file input of reads aligned to a common reference. That common reference is an assembly of pooled reads from all your experiments.

[1]: https://bitbucket.org/berkeleylab/metabat/overview
[2]: https://concoct.readthedocs.io/en/latest/
Question: Hello!
I obtained a list of unmapped reads IDs from my BAM file and I want to remap only the unmapped reads with other parameters.
How can I extract the subset of unmapped reads from my original fastq file?
Thank you in advance,
Luke

Answer: I also wrote a program for this purpose, distributed with [BBMap][1]. Usage:

filterbyname.sh in=reads.fq out=filtered.fq names=names.txt include=t

The `include` flag will toggle between including or excluding the names in `names.txt` (which can, alternately, be another fastq or fasta file). This also supports paired input/output, and names being substrings or superstrings of read IDs.

[1]: https://sourceforge.net/projects/bbmap/
Question: What are some good resources to learn shell script for NGS pipeline development?

How much shell script should one know to develop an intermediate level pipeline for NGS data analysis?

Can someone suggest some good resources, tutorials?

Answer: What is an "intermediate level pipeline"? What is your target audience? Release the pipeline into the wild? Internal lab use? Personal use? Anyway, to learn shell scripting for NGS pipeline development, you must learn shell scripting, so look at the "[Bash Guide for Beginners][1]" and "[Advanced Bash-Scripting Guide][2]".

With a very basic understanding of bash scripting you may easily put together a simple pipeline which will, for example, clean your reads, assemble a genome, map the reads / additional reads into assembled genome, and annotate assembled genome. In fact, I wrote such simple pipeline - it is really crude, no error checking, no optimizations, no whatever, but I feed fastq files and some hours later get a draft genome and its annotation.

[1]: http://www.tldp.org/LDP/Bash-Beginners-Guide/html/Bash-Beginners-Guide.html
[2]: http://www.tldp.org/LDP/abs/html/
Question: I'd want to combine 3 pheatmap into 1 figure. Here is the test code. it seems `mfrow` parameter does not work for it.

```
test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")

par(mfrow=c(1,2))
pheatmap(test)
pheatmap(test, kmeans_k = 2)
pheatmap(test, kmeans_k = 3)
```

Answer: If it's just for the one figure, I'd do it manually in Inkscape. However, for high throughput figure production, I would export individual images from R and script the final figure assembly from the components, for example using ImageMagick. You could also write an SVG file directly (for example with the perl module SVG) or by scripting Inkscape.
Question: Is there a way to convert a `.vcf.gz.tbi ` file back to a `.vcf.gz` or `.vcf` file?

My end goal is to convert these files into plink bed/bim/fam files, but plink will not accept the .vcf.gz.tbi file as is. Any help is appreciated!

Answer: No. An index (that is what tbi is) is just a table of content. You cannot reconstruct every page of a book with just the ToC at hand.
Question: I would like to parse a BAM file in parallel using pysam and multiple_iterators

Here is my code

import pysam
import sys
from multiprocessing import Pool
import time
def countReads(chrom,Bam):
count=0
#Itr = Bam.fetch(str(chrom),multiple_iterators=False)
Itr = Bam.fetch(str(chrom),multiple_iterators=True)
for Aln in Itr: count+=1

if __name__ == '__main__':
start = time.time()
chroms=[x+1 for x in range(22)]
cpu=6
BAM = sys.argv[1]
bamfh = pysam.AlignmentFile(BAM)
pool = Pool(processes=cpu)
for x in range(len(chroms)):
pool.apply_async(countReads,(chroms[x],bamfh,))
#countReads(chroms[x],bamfh)
pool.close()
pool.join()
end = time.time()
print(end - start)

I get this error when I run it.

TypeError: _open() takes at least 1 positional argument (0 given)

But it spits out a whole bunch of errors. Can anyone help me to use multiprocessing to read a BAM file in parallel using pysam?

Thanks

Answer: fixed it. I was going off some online blog that was wrong

import pysam
import sys
from multiprocessing import Pool
import time
def countReads(chrom,BAM):
count=0
# here's the fix
bam = pysam.AlignmentFile(BAM,'rb')

Itr = bam.fetch(str(chrom),multiple_iterators=True)
for Aln in Itr: count+=1

if __name__ == '__main__':
start = time.time()
chroms=[x+1 for x in range(22)]
cpu=6
BAM = sys.argv[1]
pool = Pool(processes=cpu)
for x in range(len(chroms)):
pool.apply_async(countReads,(chroms[x],BAM,))
#countReads(chroms[x],bamfh)
pool.close()
pool.join()
end = time.time()
print(end - start)
Question: One of our project used to query OMIM data as XML through NCBI's efetch utility, as described here for example:

What is the best way to interact programmatically with OMIM?

However, it seems the service has stopped functioning a few months ago. It now simply returns the following error:

Database: omim - is not supported

I can find no mention of an update to the API on NCBI's website or anywhere else.
At the same time, the pages accessible directly on OMIM's website offer no link to structured data (XML or otherwise) and the downloadable file, while using some specific format to delimit fields, is still far from the flexibility of the former XML files (for example, it is impossible to retrieve metadata for each reference).

Is there currently any way to regain access to OMIM data in a structured, parsable format (XML...)?

Answer: I've been waiting for this since July 2011.
https://twitter.com/OmimOrg/status/90770519841980416

@yokofakun[?] @mrrizkalla and we have an web API in the wings, not
quite ready to open it up yet . omim.org (@OmimOrg)
Question: Can the chromosome names and lengths to which a bam file was aligned be determined solely from its bai file?

Answer: No. The BAI format is described in §5.2 of [SAMv1.pdf](http://samtools.github.io/hts-specs/SAMv1.pdf) and does not contain the chromosome names or lengths. Instead it merely identifies reference sequences by their index, 0 <= n < n_ref, in the corresponding BAM file's (binary) header.

The CSI format operates similarly.

The [Tabix format](https://samtools.github.io/hts-specs/tabix.pdf) OTOH does contain the names of reference sequences (but not their lengths).
Question: When I searched the NCBI Taxonomy db for id=693 it has retrieved the result but the tax_id in the result is different by the ID which I searched
http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=603&retmode=xml

and I just find out this tax_id (603) is related to gi=156763600 on gi_taxid_prot.dmp.gz and If you search it on ncbi and press Taxonomy from the right part of the page it comes by error on the next page.

Have you any explanation for these!?

Answer: I assume you mistyped the first taxonomy ID and meant 603 instead. When you are using NCBI's web-search form, then a search for the ID 603 will redirect you to Salmonella enterica subsp. arizonae, which coincides with the XML query result you linked to. Said Salmonella has taxonomy ID 59203 just as in the XML result.

The taxonomy IDs 603 and 59203 appear to be denoting the same database entry, which is reflected in the XML at the bottom with:

[AkaTaxIds]
[TaxId]603[/TaxId]
[/AkaTaxIds]

I think you can treat taxonomy IDs 603 and 59203 as synonyms.

Hope that helps.
Question: I got several SAM/BAM files, and I am interested in 5'ends of the mapped reads. Is there any tools or scripts to count how many 5'ends are mapped at a specific genomic position?

N.B. I am not try to count the total reads number mapped to the specific genomic position.

For example, in the following window, I get 5 reads, 3 of them have a 5'ends at 14488, the other 2 have a 5'ends at 14487.

And I want produce a table like this.

.

Answer: bedtools genomecov -d -5 -ibam input.bam -g genome.fa > output.txt

More options see bedtools documents at http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html
Question: New version of BCFTools (using call instead of view) doesn't produce VCF files with GQs in the FORMAT field.

Does anyone know how to get it in the new version (1.1)?

Answer: When you run bcftools call specify -f GQ to have it added to the output

bcftools call -vmO z -f GQ
Question: Hi,

I've got bisulfite-sequencing data for two differentiation stages. The raw data was mapped using Bismark. For each CpG site, the methylation ratio was marked as "A/B" (A methylated reads vs B unmethylated reads for this site).

Now I want to compare the overall methylation of a certain region. Assuming this region contains 3 CpG sites, the methylated/unmethylated ratio of each site in the two stages are 45/60, 4/10, 5/15 and 14/65, 3/9, 7/20 respectively.

I think there are two ways to calculate the overall methylation rate of the region:

**Method 1: total amount of methylated reads / total number of reads within the region**

methylation rate of stage 1 = (45+4+5)/(45+60+4+10+5+15) = **54/139** = 0.39

methylation rate of stage 2 = (14+3+7)/(14+65+3+9+7+20) = **24/118** = 0.20

**Method 2: average the methylation rate of each CpG site within the region**

methylation rate of stage 1 = (45/(45+60) + 4/(4+10) + 5/(5+15))/3 = 0.32

methylation rate of stage 2 = (14/(14+65) + 3/(3+9) + 7/(7+20))/3 = 0.22

For **method 1**, I can then calculate the confidence inteval and the P value to know whether **54 out of 139** is significantly different from **24 out of 118**.

However, the methylation rate calculated by **method 1** is biased toward the site with higher reads coverage.

**Method 2** seems to be more robust to indicate the actual methylation rate of the region.

But I don't know which statistics test should I use for method 2 to know the significance.

Please help.

Thanks in advance.

Answer: From context I gather that you lack biological replicates, which is unfortunate but still rather common. As you mentioned, "method 1", which is basically summing across CpGs and performing a Fisher's test is less than ideal, since it will miss interesting cases (using your earlier nomenclature, think of 50/0, 0/50 in sample1 and 25/25, 25/25 in sample 2, they're very different but this method would miss that).

Method 2 lends itself to a weighted paired t-test, which you can do in R. This method is, of course, also not ideal, since you give a weight to each CpG, but if coverage is drastically different between the two samples at a CpG then you start running into problems (there are ways around that, but they'd be moderately annoying).

What would likely be a nice method in this case would be to perform a Fisher's test on the individual CpGs between the samples and then spatially correct the p-values (the original paper from Benjamini & Heller is here but you can find an implementation in the BiSeq bioconductor package, which you could also just directly use to make your life easier since it's actually targeted at RRBS or other targeted BS-seq datasets). You could also use the SLIM model for p-value adjustment, such as is used in methylKit. Somewhat similarly, you could also first use a smoothing method (see BSseq as an example) and then do a paired t-test, though I wouldn't recomment that for your use case.

There are some additional possibilities, but I'd recommend just using BiSeq and being done with it.
Question: I have a list of genomic ranges mapped to hg19 . my data is in matrix format lets call it `ranges` which has 600,000 rows and 4 clumns
here is few row of my data

head(ranges)
chr start end strand
[1,] "chr1" "10025" "10525" "."
[2,] "chr1" "13252" "13752" "."
[3,] "chr1" "16019" "16519" "."
[4,] "chr1" "96376" "96876" "."
[5,] "chr1" "115440" "115940" "."
[6,] "chr1" "235393" "235893" "."
Is there a function that gets sequences and calculates GC content for each row( each range)
I would prefer that output be in a vector format
I would really appreciate your help

Answer: thank you all,
these solutions are great.
I also found a function that calculates GC content using BSgenome.Hsapiens.UCSC.hg19.

require(BSgenome.Hsapiens.UCSC.hg19)
require(BSgenome)
library(Repitools)
Granges<- makeGRangesFromDataFrame(data.frame(ranges))
gc<- gcContentCalc(Granges , organism=Hsapiens, verbose=TRUE)
Question: I have many files like this

ENSG00000000003.13	366
ENSG00000000005.5	26
ENSG00000000419.11	1905
ENSG00000000457.12	775
ENSG00000000460.15	377
ENSG00000000938.11	316
ENSG00000000971.14	272
ENSG00000001036.12	1726
ENSG00000001084.9	2479
ENSG00000001167.13	1166
ENSG00000001460.16	38
ENSG00000001461.15	298

They are all htseq-count.

I have two questions

1- how can I convert the ID to gene names? do you have any solution available (I can use any language) ?
2- what would you use for normalizing them?

Answer: To answer your first question :

Create a TSV file with your input => list.tsv

Remove ".13" etc from your ENSG (list.tsv) => list_ensembl_gid.tsv

awk -F '\t' -v OFS='\t' '{sub(/\.[0-9]*/, "", $1)} 1' list.tsv > list_ensembl_gid.tsv;

Then use R & Biomart to retrieve corresponding genes names

R
data = read.table("list_ensembl_gid.tsv")
colnames(data)[1] <- "ensembl_gene_id"
colnames(data)[2] <- "counts"

library('biomaRt')

hsapiens = useMart("ensembl",
dataset="hsapiens_gene_ensembl")

hsapiens_infos <- getBM(attributes=c('ensembl_gene_id',
'external_gene_name'),
mart = hsapiens)

merge_infos <- merge(x = data,
y = hsapiens_infos,
by = "ensembl_gene_id",
all.x = TRUE)

head(merge_infos)

> ensembl_gene_id counts external_gene_name

> ENSG00000000003 366 TSPAN6

> ENSG00000000005 26 TNMD

> ENSG00000000419 1905 DPM1

> ENSG00000000457 775 SCYL3

> ENSG00000000460 377 C1orf112

> ENSG00000000938 316 FGR
Question: This is a very basic question which I have not been able to understand even after searching on the web for significant amount of time. Can someone please explain me as if to a high-school student about these terms:

1. contig
2. Singleton read; and
3. Discordant mate pairs

Answer: Contigs are long continuous stretch of sequence results from overlapping short reads.

When you do a paired-end sequence, you sequence a DNA fragment ( lets say of length 400bp ) from both the ends, resulting R1 and R2 reads. When you map them back to genome, if the both reads map with in 400bp (as they are generated from 400bp fragment) and expected strand orientation, they are referred as concordant pairs. This is primarily in DNA-Seq. I don't think there is a discordant mate pair, it is called discordant pair. Mate-Pair sequencing is a different thing.

When one of the mates (either R1 or R2) is not present (lost due to pre-processing or in mapping step or post-processing of data), the other mate becomes **single (also sometimes referred as orphan reads) called as singletons**.
Question: Does anyone know how to solve this error? I have never seen this error before and my current script has been working fine up until this point.

> t2g <- biomaRt::getBM(
+ attributes = c("ensembl_transcript_id", "transcript_version",
+ "ensembl_gene_id", "external_gene_name", "description",
+ "transcript_biotype"), mart = mart)
No encoding supplied: defaulting to UTF-8.
Error in biomaRt::getBM(attributes = c("ensembl_transcript_id", "transcript_version", :
The query to the BioMart webservice returned an invalid result: biomaRt expected
a character string of length 1. Please report this to the mailing list.

Answer: For me the Ensembl mirror sites are down. It looks like you're US based, so by default you'll get redirected to [uswest.ensembl.org](http://uswest.ensembl.org/biomart/martview?redirect=no) or [useast.ensembl.org](http://useast.ensembl.org/biomart/martview?redirect=no), neither of which are working for me at the moment.

The main site is working so you want to use [www.ensembl.org][1] as your host. However you also have to supply the argument `ensemblRedirect = FALSE`, otherwise their internal redirection will simply send you back to your local site. This works for me, even if I'm using a server based in Texas:

mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "hsapiens_gene_ensembl",
host = 'www.ensembl.org',
ensemblRedirect = FALSE)

After some discussion with the Ensembl BioMart team the plan going forward is to remove the redirection entirely when using **biomaRt**, so this should no longer be an issue, and you'll go to whatever address you've provided. I'll try to update here when that's been done.

----------

The 'No Encoding Supplied' message is a red herring relating to the content being returned by the server - you see it whether there's a problem or not. I've already addressed this in the developmental version of **biomaRt**, and will update the release version too. It's not very helpful to end users!

[1]: http://www.ensembl.org?redirect=no
Question: Hi,

I've been using the blastn (version 2.2.28+) standalone tool against a custom formatted genome via:

blastn -db BLASTDB -word_size 7 -query input.fa -out filename -perc_identity 100 -outfmt 6 -max_target_seqs 2

To discard non-perfect hits and show only the 2 top hits.

The output file has a great format however is there a way to add an extra column that contains the actual target-seq (sequence of the matched hit)? Such that the fields are:

query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score, sequence

Thanks!

TJC

Answer: Run `blastn -help` then look for the field called `outfmt`

*** Formatting options
-outfmt
alignment view options:
0 = pairwise,
1 = query-anchored showing identities,
2 = query-anchored no identities,
3 = flat query-anchored, show identities,
4 = flat query-anchored, no identities,
5 = XML Blast output,
6 = tabular,
7 = tabular with comment lines,
8 = Text ASN.1,
9 = Binary ASN.1,
10 = Comma-separated values,
11 = BLAST archive format (ASN.1)

Options 6, 7, and 10 can be additionally configured to produce
a custom format specified by space delimited format specifiers.
The supported format specifiers are:
qseqid means Query Seq-id
qgi means Query GI
qacc means Query accesion
qaccver means Query accesion.version
qlen means Query sequence length
sseqid means Subject Seq-id
sallseqid means All subject Seq-id(s), separated by a ';'
sgi means Subject GI
sallgi means All subject GIs
sacc means Subject accession
saccver means Subject accession.version
sallacc means All subject accessions
slen means Subject sequence length
qstart means Start of alignment in query
qend means End of alignment in query
sstart means Start of alignment in subject
send means End of alignment in subject
qseq means Aligned part of query sequence
sseq means Aligned part of subject sequence
evalue means Expect value
bitscore means Bit score
score means Raw score
length means Alignment length
pident means Percentage of identical matches
nident means Number of identical matches
mismatch means Number of mismatches
positive means Number of positive-scoring matches
gapopen means Number of gap openings
gaps means Total number of gaps
ppos means Percentage of positive-scoring matches
frames means Query and subject frames separated by a '/'
qframe means Query frame
sframe means Subject frame
btop means Blast traceback operations (BTOP)
staxids means Subject Taxonomy ID(s), separated by a ';'
sscinames means Subject Scientific Name(s), separated by a ';'
scomnames means Subject Common Name(s), separated by a ';'
sblastnames means Subject Blast Name(s), separated by a ';'
(in alphabetical order)
sskingdoms means Subject Super Kingdom(s), separated by a ';'
(in alphabetical order)
stitle means Subject Title
salltitles means All Subject Title(s), separated by a '<>'
sstrand means Subject Strand
qcovs means Query Coverage Per Subject
qcovhsp means Query Coverage Per HSP