cannin/biostars_qa huggingface dataset:
Question: Hi all,

For doing GWAS and eQTL analyses typically many covariates are included. Usually known covariates like gender, age etc. and also unknown covariates like PCA dimensions. Problem is some of these covariates might be not independent and have a relatively high correlation. Does anybody knows a what a common cut of value is for correlation or what to do with not independent covariates like age and PCA dimensions?

Answer: I prefer to avoid using highly correlated covariates as they can often be redundant.

My approach is to use the following procedure:

1. Select the (non-PCA) covariates that you are certainly including in the analysis.
2. Perform PCA analysis on the data with the covariates selected in step 1 regressed out.
3. Note that top PCs from step 2 are guarantied to be orthogonal to the covariates selected in step 1.
4. Run eQTL analysis with covariates from step 1 and a few top PCs from step 2.
Question: Hi,

to evaluate genome assemblers I would like to simulate PacBio HiFi reads. However, there seems no read simulator available that has a HiFi mode (at least I did not find anything via Google. If someone can point me to something, it would be amazing.).

So I would just write my own simple read simulator for starters, and I was thinking about just adding random indels and subs with a 1% rate. This would roughly match the minimum Q20 that PacBio defines their HiFi reads as.

Is this a viable approach, or would I be missing some important characteristic?

Answer: You could try something like I did here: https://github.com/chhylp123/hifiasm/issues/33
Question: Hi, I have output from snp-dists (https://github.com/tseemann/snp-dists) in molten format, e.g.:

seq1 seq2 1
seq1 seq3 2
seq2 seq1 1
seq2 seq3 3
seq3 seq1 2
seq3 seq2 3

The third column gives the number of SNPs between the pair of sequences given in columns 1 and 2. As you can see, these values are duplicated, as it shows both the combination seq1 seq2 and seq2 seq1. How can I (in R or bash preferably) remove the duplicate values?

Answer: Using **awk**:

$ awk '!(seen[$1,$2]++ || seen[$2,$1]++)' test.txt
seq1 seq2 1
seq1 seq3 2
seq2 seq3 3

Using **R**:

# example data
x <- read.table(text = "seq1 seq2 1
seq1 seq3 2
seq2 seq1 1
seq2 seq3 3
seq3 seq1 2
seq3 seq2 3")

# sort column values, then get unique
unique(data.frame(c1 = pmin(x$V1, x$V2), c2 = pmax(x$V1, x$V2), value = x$V3))
# c1 c2 value
# 1 seq1 seq2 1
# 2 seq1 seq3 2
# 4 seq2 seq3 3

Using **R** again, a bit simpler and scales better when we have more than 2 columns, ([Related StackOverflow post](https://stackoverflow.com/q/9028369/680068)):

x[ !duplicated(apply(x[, 1:2], 1, sort), MARGIN = 2), ]
Question: Hi,

I am trying to do differential expression based on RNA-seq data. I only have one variable, which is the timepoint. I want to analyze gene expression relative to the 0 h time point. I already set up the contrasts. Now, I am trying to set up the design matrix. What is the difference between

design <- model.matrix(~timepoint)

and

design <- model.matrix(~0 + timepoint).

There is a big difference in the results.

Answer: The clue here is to look at the coefficients, or columns derived in the design matrix. Your non intercept model design, will create a coefficient that describes the difference between the two levels of your `timepoint` variable. Your non-intercept model will create a coefficient for each level. Here's an example using your non-intercept model:

library(tidyverse)
library(limma)
# Make a Phenotype Table
foo <- data.frame(SampleType = paste0(rep(c("A","B"), each = 3)),
Reps = rep(1:3, 2)) %>%
mutate(Names = paste0(SampleType,Reps))
# Make some fake gene expression data
genes.foo <- matrix(rnorm(6*100),ncol = 6) %>%
`rownames<-`(paste0("Gene_",1:100)) %>%
`colnames<-`(foo$Names)

#non-intercept model
model.foo <- model.matrix(~0 + SampleType, data = foo) %>%
`colnames<-`(gsub("SampleType","",colnames(.)))

model.foo
A B
1 1 0
2 1 0
3 1 0
4 0 1
5 0 1
6 0 1

conts.foo <- c("A_Vs_B" = "A-B")
contmap.foo <- makeContrasts(contrasts = conts.foo, levels = colnames(model.foo)) %>%
`colnames<-`(names(conts.foo))
fit.foo <- lmFit(genes.foo, model.foo) %>%
contrasts.fit(contmap.foo) %>%
eBayes

> fit.foo$coefficients %>% head
Contrasts
A_Vs_B
Gene_1 -0.6722254
Gene_2 -1.1571439
Gene_3 -0.1021150
Gene_4 0.8687436
Gene_5 -0.9600460
Gene_6 1.3139228

topTable(fit.foo, coef = "A_Vs_B", number = Inf, p.value = 0.05, lfc = log2(1.5))

Alternatively, here's your intercept model

model.foo <- model.matrix(~SampleType, data = foo)

fit.foo <- lmFit(genes.foo, model.foo) %>%
eBayes

> model.foo
(Intercept) SampleTypeB
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1

# Here SampleTypeB is actually the difference relative to the first level (SampleTypeA)
> fit.foo$coefficients %>% head
(Intercept) SampleTypeB
Gene_1 -0.2920748 0.6722254
Gene_2 -1.0356549 1.1571439
Gene_3 0.2477702 0.1021150
Gene_4 0.5016585 -0.8687436
Gene_5 0.3153554 0.9600460
Gene_6 0.2836491 -1.3139228

topTable(fit.foo, coef = "SampleTypeB", number = Inf, p.value = 0.05, lfc = log2(1.5))

edit: Changed the data frame initialisation.
Question: Hello,

Does cellranger count require index fastqs? The doc specifies the output directory of cellranger mkfastq (or csv options), however I have worked with the R1/R2 fastqs generated from another lab without issue. Thanks.

-Todd

Answer: You have likely [seen this page][1] on 10x support site. Sample index files are optional but you do need to have your datafiles named in a specific way (`bcl2fastq`). See section on `My FASTQs are not named like any of the above examples` at end of the page linked.

[1]: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input
Question: I did a reciprocal BLAST search to identify homologues genes (orthologs and paralogs) from 2 different species (A, B) in the genome of species C. What I'm trying to find out is whether e-values of homologs from species A are significantly higher then from species B, and whether A is more likely to have orthologes/paralogs in C.

Is that even possible or does it make any sense? How do I calculate that? Another problem is that I also have E-values of 0 (< 2.225074E-308 I believe) and I don't know how to deal with them in a statistical evaluation (obviously the best hits).

I am rather new to bioinformatics ... and please excuse my english.

Thanks

Answer: As Jean-Karim pointed, e-values are not comparable between BLAST searches of two different databases. Bitscore won't help you either, because this score also depends on the size of the database.

**Solution 1 (more accurate):**

Instead of playing with BLAST scores by yourself, just use the [InParanoid][1] software, which finds orthologs and paralogs between two FASTA files (A and C). InParanoid is based on Reciprocal Best BLAST Hit, but it does a lot of things for you. As a result, you get a list of all orthologs and paralogs and, most importantly, confidence value for each ortologous/paralogous pair. In this way, you can apply t-test (if values are normally distributed) or Mann-Whitney test (if they are not).

**Solution 2 (primitive, questionable):**

If you don't want to use InParanoid, just simple RBH then you need to tell BLAST that the search space is equal among BLAST of the same query to multiple databases. Then, both, E-values and bitscores will be comparable among all results because the search space will be consistent. You do this by simply using `-dbsize` flag when running BLAST locally. In my opinion, you can use the average of the sequences in A, B and C. For example, if A has 4000 sequences, B has 3000 and C has 5000 sequences, I would set -dbsize flag to 4000.

For example:

blastp -query A.fasta -db C.fasta -dbsize 4000

[1]: http://software.sbc.su.se/cgi-bin/request.cgi?project=inparanoid
Question: Hi !

I have many .bam that I want to get their .bai using samtools in the terminal.

I tried the following command :

samtools index *.bam

However, I did not get any .bai file.

Regards

Answer: using GNU parallel:

parallel samtools index ::: *.bam
Question: Hi,

running this

> ensembl = useMart(biomart
> ="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl", host="may2017.archive.ensembl.org")

or this:

> ensembl = useMart(biomart
> ="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl", host="http://dec2017.archive.ensembl.org")

used to work in the past. but now it gives the error:

> Request to BioMart web service failed. Verify if you are still
> connected to the internet. Alternatively the BioMart web service is
> temporarily down. Check http://www.biomart.org and verify if this
> website is available. Error: XML content does not seem to be XML:

The site is up, the internet connection also.
In the site it says this:

> If you are a user of biomaRt (a part of the Bioconductor library)
> change the host from 'www.biomart.org' to 'www.ensembl.org'

But in the commands I used to run there's no 'www.biomart.org' to replace...

Do you know what is the correct command now?

I also tried:

> ensembl =
> useMart(biomart="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl",
> host="www.ensembl.org")

and didn't work

Answer: As Emily has said, parts of Ensembl are down right now. This includes the archive sites and the BioMart interface at the main site. You can check this just by visiting in a browser (http://www.ensembl.org/biomart/martview?redirect=no) - at the moment you get an error page.

However you can use one of the mirror sites e.g. http://uswest.ensembl.org/biomart/martview?redirect=no

To do this using the **biomaRt** package you would do something like:

ensembl = useMart(biomart="ENSEMBL_MART_ENSEMBL",
dataset="hsapiens_gene_ensembl",
host="uswest.ensembl.org",
ensemblRedirect = FALSE)

I would also suggest updating your version of **biomaRt** (and maybe also R). The error message about visiting www.biomart.org is misleading and is no longer produced in more recent versions of **biomaRt**. You now get a more appropriate URL to try, that would have taken you to the error page at Ensembl.
Question: Hi,

I aligned a few samples using STAR to the genome provided in the Illumina iGenomes UCSC hg19 bundle ([here][1]) -- I used the provided gene feature (gtf2) file as is. Now, my motive is to calculate the gene and isoform expression levels using bedtools multicov (at the same time).

Use of the gtf2 file produces a file containing read counts per exon. I wish to compute gene and isoform read counts too, so I converted the gtf2 file to a gff3 file using using gtf2gff3 script from SO/GAL ([here][2]).

My first question is: Is it OK if the alignment is performed with gtf2 file but counted for reads using the gff3 file, keeping in mind that the gff3 file was converted from the gtf2 file?

My second question follows

I have read both these resources ([here][3] and [here][4]) but do not understand the differences between:

- exon vs CDS
- transcript vs mRNA

I know that with the process I described, it is possible to retrieve gene read count by selecting only the lines where feature=gene from the bedtools multicov output. What must I do for isoforms? I am confused by the semantics.

Thanks ahead of time and let me know if my post was not clear enough.

[1]: http://cufflinks.cbcb.umd.edu/igenomes.html
[2]: http://www.sequenceontology.org/software/GAL.html
[3]: http://mblab.wustl.edu/GTF2.html
[4]: http://www.sequenceontology.org/gff3.shtml

Answer: Regarding mRNA vs. transcript, since Istvan already mentioned exons and coding sequences (CDS), mRNAs are a subset of transcripts that encode proteins. There are many other types of transcripts (tRNAs, miRNAs, lncRNAs, etc.) that are not protein-coding. Depending on the type of library that was sequenced, you may expect to find mostly mRNAs or mostly other forms.
Question: Hi,

I have one .vcf file of whole genome sequencing of tumour Vs normal samples of 21 patients.

I need a data from like this as input for a tool for finding driver genes

> head(mutations)
sampleID chr pos ref mut
1 Sample_1 1 871244 G C
2 Sample_1 1 6648841 C G
3 Sample_1 1 17557072 G A
4 Sample_1 1 22838492 G C
5 Sample_1 1 27097733 G A
6 Sample_1 1 27333206 G A

In separated .vcf files for each patient I have start, end, chromosome, ref, and variant allele. However I am sure how to get such data frame from this big vcf

Any help please?

Thank you

Answer: Using *bcftools*:

bcftools query -f '[%SAMPLE %CHROM %POS %REF %ALT %GT\n]' myFile.vcf > myFileLong.txt
Question: I've got a list of rs SNPs I'd like to enter into the David functional annotation tool. Since it does not support rsids I need to get the refseq gene id (or something similar) of the genes these SNPs overlap with first.

Is there a simple way of getting a gene ID for a SNP? Solution in BioPython or R is fine.

Since I have 22k snps I need an automatic way. And if a text file that maps these values exist, a link to it would be enough.

Ps. preferably looking for a solution that does not use the position of the SNPs; I would be able to solve the problem this way myself.

Answer: It looks like DAVID can convert from Ensembl Gene IDs, so you can get to those from rs IDs using [R/biomaRt][1], like this:

```
library(biomaRt)
mart.snp <- useMart("snp", "hsapiens_snp")

getENSG <- function(rs = "rs3043732", mart = mart.snp) {
results <- getBM(attributes = c("refsnp_id", "ensembl_gene_stable_id"),
filters = "snp_filter", values = rs, mart = mart)
return(results)
}

# default parameters
getENSG()
refsnp_id ensembl_gene_stable_id
1 rs3043732 ENSG00000175445

# or supply rs ID
getENSG(rs = "rs224550")
refsnp_id ensembl_gene_stable_id
1 rs224550 ENSG00000262304
2 rs224550 ENSG00000196689
```

[1]: http://www.bioconductor.org/packages/release/bioc/html/biomaRt.html
Question: Hello,

I would like you to suggest me a literature-based interactomics tool. Thanks in advance.

Answer: iHop is a mature system that uses gene name recognition text-mining on MEDLINE abstracts and allows you to build a gene/protein network model by clicking and adding sentences with multiple gene names. From the iHop help pages:

Nodes in this graph represent genes.
Edges correspond to sentences that
associate two genes with each other...The model graph
makes it possible to analyze your
collection of sentences in an
interactive manner and to get familiar
with the newly aquainted knowledge. By
clicking on a gene, all its synonyms
will be highlighted in red and the
corresponding sentences will be ranked
to the top. By clicking on an edge,
those sentences containing the
connected genes are shown at the top.
Furthermore, a line separates these
highlighted sentences from all others...By clicking on a
gene in a sentence, the corresponding
information page for the gene will be
opened and additional information can
be added to the Gene Model.

Another option that has been around for some time is PubGene:

PubGene mines the abstract texts of 25 Million PubMed articles for co-citation of multiple genes or proteins and displays them as "Literature Networks", where nodes represent each gene or protein and the connecting lines represents the number of articles, in which each gene or protein pair is co-cited....When a gene or protein is studied, there is a good chance its name (or a synonym for that name) will appear in articles together with other gene or protein names. One can visualize how most genes that have been studied will be connected either directly or indirectly to each other in a Literature Network. Connections in the literature are a strong indicator of biological interaction.
Question: I have fasta file namely `119XCA.fasta` as shown below,

>cellulase
ATGCTA
>gyrase
TGATGCT
>16s
TAGTATG
I need to remove all the fasta headers, keep the sequences one by one and need to write file name as a fasta header. The expected outcome is shown below,

>119XCA
ATGCTA
TGATGCT
TAGTATG
I have used the following script `sed '/^>/d' foo.fa > out.fa` which remove the fasta headers but, i do not know how to manage to write file name as a header. Therefore, please help me to do the same.

Answer: try this:

$ cat test.fa
>cellulase
ATGCTA
>gyrase
TGATGCT
>16s
TAGTATG

$ awk 'BEGIN {print ">"ARGV[1]};!/^>/{print}' test.fa

>test.fa
ATGCTA
TGATGCT
TAGTATG

$ cat <(echo ">"$basename test.fa) <(grep -v ">" test.fa) (note:extra space in header)
> test.fa
ATGCTA
TGATGCT
TAGTATG
Question: Project PRJEB99111 has 147 samples. I want to download the metadata (age, sex, disease status, etc) of each sample, not fastq. The only way I can download the metadata is by downloading the xml file of each sample accession one by one - is there a way to bulk download all 147 metadata files? I can work with xml files if I have to.

You can view the metadata for a specific sample accession by clicking on the"attributes" tab. Here is an example for one sample: [https://www.ebi.ac.uk/ena/data/view/SAMEA104228123][1]

[1]: https://www.ebi.ac.uk/ena/data/view/SAMEA104228123

Answer: Using NCBI eUtils: `esearch -db bioproject -query "PRJEB99111" | elink -target biosample | efetch -format docsum | xtract -pattern DocumentSummary -block Attribute -element Attribute`

produces (only a sample below)

2017-08-28	2017-08-26	ERS1887138	female	44	years	UBERON:feces	UBERON:feceUBERON:feces	UCSF	2013-09-25	Missing: Not provided	RRMS	MS	TRUE	124	urban biome	human-associated habitat	feces	human-gut	No	USA:CA:San Francisco	Missing: Not provided	Missing: Not provided	Homo sapiens	111	9606	Missing: Not provided	mimarks-survey	37.76	adult	-122.46	UCSF	FALSE	1	stool	1	Missing: Not provided	1	1_a	No	Gut dysbiosis in patients with multiple sclerosis is characterized by bacteria that regulate T lymphocyte differentiation in vitro	No_Treatment	Off	Missing: Not provided	Missing: Not provided	dry	1990
2017-08-28	2017-08-26	ERS1887137	male	61	years	UBERON:feces	UBERON:feceUBERON:feces	UCSF	Missing: Not provided	Missing: Not provided	RRMS	MS	TRUE	124urban biome	human-associated habitat	feces	human-gut	No	USA:CA:San FranciscMissing: Not provided	Missing: Not provided	Homo sapiens	62	9606	Missing: Not provided	mimarks-survey	37.76	adult	-122.46	UCSF	FALSE	1	stool	2	Missing: Not provided	1	1_a	No	Gut dysbiosis in patients with multiple sclerosis is characterized by bacteria that regulate T lymphocyte differentiation in vitro	No_Treatment	Off	Missing: Not provided	Missing: Not provided	dry	1984
Question: Hi,

I am new to Bioinformatics.I am trying to use GATK tool for finding SNP and indels but the problem is that the documentation seems to be complex for a beginner like me with so many tools to start with including IndelRealigner,MuTect,Haplotypecaller..The more I read the more confused I am getting.I also tried using Indel religner which is the first step in the process of identification of indels(correct me if wrong) but I got an error saying that you don't have the intervals file in the targetIntervals in the following command.

java -jar GenomeAnalysisTK.jar \
-T IndelRealigner \
-R reference.fasta \
-I input.bam \
-known indels.vcf \
-targetIntervals intervalListFromRTC.intervals \
-o realignedBam.bam

The error message is:Could not read file /home/aditya/Bioinfotools/gatk/intervalListFromRTC.intervals because The interval file does not exist.
I don't know where to get the intervals file from.

But please guide me overall in understanding GATk.

Answer: There is a step by step guide: GATK best practices
https://www.broadinstitute.org/gatk/guide/best-practices.php
Question: Hi, everybody. I have FASTQ headers of the form

@FCC3KD2ACXX:6:1101:1545:2184#ATCACGATC/1

The "ATCACGATC" portion of these older-style headers is supposed to be the "index sequence", or the molecular barcode of a multiplexed sample, according to the Wikipedia article. But I know what the barcodes are, and that isn't one of them. All the barcodes for this project are 6-8bp, not 9, and there are "only" about 300 legitimate barcodes in the lane.

Overall, there are over 255,000 different individual ones of these index sequences on various different reads, out of about 178M reads in the file. Some of them even contain Ns, but they're all exactly 9bp long. And this particular sequence (ATCACGATC) is by far the most prevalent -- it's on 90% of the reads, so it can't be succeeding at separating anything out very specifically.

I'm coming in late to a project, and all I have to go on is the FASTQ files, the list of barcodes, and some wet-lab protocol docs that I'm not particularly qualified to interpret. Any idea what this odd extraneous-looking sequence is? If so, thanks in advance!

Answer: Your sequence string 'ATCACGATC' is a perfect match to the TruSeq adapter index 1. The length (9bp) is atypical but not unheard-of. The length of the index read is often index-length-plus-one (i.e., seven cycles for a 6mer), so it appears that the sequencer was programmed for nine index cycles (8mer+1). The last nucleotide is typically trimmed from the data by default (--use-bases-mask = I8n), but the full length can be specified by the user (I9).

As @igor indicates, the amount of data suggests that the file was not demultiplexed. However, given that 90% of your data are a perfect match, it's also unlikely to contain multiple libraries. Best guess is that some other lanes of this flow cell did contain 8mer indices, and CASAVA doesn't allow lane-specific BCL-to-FASTQ conversion (so the data in all lanes have identical formats).
Question: Hi all,

in relation to a mail from January this yearin the r-help community, I followed Simon's advice to do my analyses in DESeq2 instead of DESeq.

I am working on an RNASeq from c. elegans. I have mapped the data with the ensembl genome build WBcel215. I have ran tophat2 to map and featureCounts to counts the reads (both with the defaults parameters).

I have two conditions, control and a knock-out with each three replica. Now I am trying to find differentially regulated genes between the two conditions using DESeq2.

This is the script I am using to read my raw count table into DESeq2:

```
featureCountTable <- read.table("featureCountTable_RawCounts.txt", sep="\t", quote=F)
colData <- data.frame(row.names=names(featureCountTable), condition = c(rep("wt",3), rep("cpb3", 3)))
cds <- DESeqDataSetFromMatrix (
countData = featureCountTable,
colData = colData,
design = ~ condition
)
fit = DESeq(cds)
res = results(fit)
```

But I am getting the same problem with DESeq2 as I have got with DESeq. When I ran the DESeq command I get a warning:
```
Warning messages:
1: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced
2: step size truncated due to divergence
```

So again I have tried to change the fitType.

fit = DESeq(cds, fitType="local")

Which than came back without any warnings.

Apparently this time both fitTypes are almost similar (at least to my inexperienced eyes.)

I add both dispersion Plots. The red line goes through the point-cloud in both cases (as Simon defined a good fit in the last communication, I wish it would have bin so easy :-) .

In the local fit type there a more outliers and the right end of the slope is going up again. I am not sure whether or not this is a good thing or not.

So, my question is - which of the two options is better?

I understand, that in general the parametric (default) option is better, but here it gives me a warning, so that something in the fit calculations is not good.

How can I understand theses plots?

Thanks for the help,

Assa

default/parametric fit

local fit

P.S.

I also tried to post it on the bioC help site, but got no responses, so I try it here.

Answer: This warning message can be ignored. It is coming from a call to R's glm() function in capturing the (dispersion ~ mean) trend. And the trend is fit iteratively until convergence, so though glm() complained at some step, it produced a final fit without error. If the parametric trend does not converge, local fit is substituted.

What version of DESeq2 are you using? I thought I had worked on more comprehensible warning reporting in this function.

I would go with the parametric to avoid the curve at the right side, although it shouldn't matter much.