cannin/biostars_qa huggingface dataset:
Question: Hi all,

I have a table that you will see below:

CHROM_POS A_Freq M.F Annotation N_Chr POP
1 CM009840.1_932 1.000000 0.000000 nongenic 20.00000 KHUZ
2 CM009840.1_1096 0.666667 0.333333 nongenic 13.33334 KHUZ
3 CM009840.1_1107 0.277778 0.277778 nongenic 5.55556 KHUZ
4 CM009840.1_1177 0.500000 0.500000 nongenic 10.00000 KHUZ
5 CM009840.1_1276 0.555556 0.444444 nongenic 11.11112 KHUZ
6 CM009840.1_1295 0.555556 0.444444 nongenic 11.11112 KHUZ
7 CM009840.1_1518 0.937500 0.062500 nongenic 18.75000 KHUZ
8 CM009840.1_1527 0.000000 0.000000 nongenic 0.00000 KHUZ
9 CM009840.1_1533 0.937500 0.062500 nongenic 18.75000 KHUZ
10 CM009840.1_1630 0.062500 0.062500 nongenic 1.25000 KHUZ

SO, I want to draw a plot like the following plot:

![enter image description here][1]

[1]: http://uupload.ir/files/4lxb_a_freq.png

What is the best idea?

Answer: Here is a start:

library(dplyr)
library(ggplot2)

# example data
set.seed(1); myData <- data.frame(
A_Freq = runif(1000),
Annotation = sample(LETTERS[1:3], 1000, replace = TRUE))

# prepare data, use "cut" make "A_Freq" groups
plotDat <- myData %>%
mutate(AlleleFrequency = cut(A_Freq, seq(0, 1, 0.25))) %>%
group_by(AlleleFrequency, Annotation) %>%
summarise(FractionOfSNPs = n()/nrow(myData) * 100)

# then plot
ggplot(plotDat,
aes(AlleleFrequency, FractionOfSNPs, group = Annotation, col = Annotation)) +
geom_line() +
scale_y_continuous(limits = c(0, 100))
Question: I am trying to download entire dataset for a bioproject using esearch and efetch from the Entrez Utilities.

My syntax is based on syntax posted by @Istvan Albert at https://www.biostars.org/p/111040/#359440, which is

> `esearch -db sra -query PRJNA40075 | efetch --format runinfo | cut
> -d
> ',' -f 1 | grep SRR | head -5 | xargs fastq-dump -X 10 --split-files`

For the BioProject PRJNA269201 I am interested in, slightly truncated syntax as shown below, creates 144 empty files as expected:

esearch -db sra -query PRJNA269201 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | xargs touch

However, when I try the full-length syntax, it behaves differently from what I expected under both scenarios 1 and 2 detailed below:

**Scenario 1**. On head-node of a cluster:

esearch -db sra -query PRJNA269201 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | head -2 | xargs fastq-dump --split-files

one file finished download, but it is 5.5G which is way larger than the 1.2GB I expected based on info at this [link][1] - is this difference because of file compression?! How can I download to a much more compressed version for both storage and downstream RNA-Seq analyses?

-rw-rw-r-- 1 aksrao aksrao 1.1G Jan 19 19:47 SRR1726554_1.fastq

-rw-rw-r-- 1 aksrao aksrao 5.5G Jan 19 19:44 SRR1726553_1.fastq

**Scenario 2**. When I try to submit this as a shell script, the STDERR stream (SLURM queue management on UBUNTU cluster) captures the following error message:

> `2019-01-20T02:28:55 fastq-dump.2.8.2 err: param empty while validating`
> `argument list - expected accession`

This same problem was reported on the original post by user @ bandanaschapagain, but it may not have been answered and resolved, hence I am posting this afresh. Could someone please help me? Thank you!
[1]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR1726553

Answer: Download the [RunInfo][1] table and use parallel to download multiple files at once.

#!/bin/bash
#change the number after -j change the number of files to be processed.
parallel --verbose -j 20 prefetch {} ::: $(cut -f5 SraRunTable.txt ) >>sra_download.log
wait
parallel --verbose -j 20 fastq-dump --split-files {} ::: $(cut -f5 SraRunTable.txt ) >>sra_dump.log
wait
exit

[1]: https://www.ncbi.nlm.nih.gov/Traces/study/?WebEnv=NCID_1_74280305_130.14.18.97_5555_1548224368_3125055013_0MetA0_S_HStore&query_key=3
Question: Hello,

I am carrying out a metagenomic study and I would like to know if there are studies listing microbial species associated with environmental contamination. If I find, for instance, crAssphage in stools, this virus is reported as a common presence in the intestinal locus. Is this virus present as a contaminant let's say on a bench as well?

Are there species that signal contamination of samples?

Answer: Yes, metagenomics is full of contaminating sequences and there is a whole literature on this going back many years, just a few of which are represented here:

- the worse the assembly, the greater the contamination problems
- contaminants are technology-specific nanopore has different contaminants to illumina
- many labs have contaminants too
- clean environments and positive / negative controls are a must

https://www.biorxiv.org/content/10.1101/2020.01.26.920173v1.full.pdf Salzberg 2020 - RefSeq and Genbank contaminants

Infamous species (mostly with illumina adapters)
Achromobacter
Turkey
various band worms
Carp Cyprinus carpio

https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-019-0678-6

Reagent microbiome? https://www.nature.com/articles/s41564-018-0202-y

Quality and Contaminant free reference genomes?

FDA-ARGOS
https://www.nature.com/articles/s41467-019-11306-6

Practically removing contaminants: https://www.molecularecologist.com/2017/01/handling-microbial-contamination-in-ngs-data/
Question: Hi everyone

I need to locate a blast local database that was created on s server long time ago, but i don't know what is it extension or how it looks if I see it through more command.

How I can find this file?

Answer: At its simplest form, blast databases have three files: for nucleotides, file.nhr, files.nin, file.nsq; for proteins, file.phr, files.pin, file.psq. Search for one of these extensions:

find . -name "*.phr"
Question: I have the following RNA-seq set-up:

- 5 samples for comparison
- matrix with TMM-normalized FPKM values

Now, TMM-normalization allows comparison across samples, so I can compare gene X between samples 1-5. Can I also compare gene X with gene Y within, say, sample 2 or does the TMM-normalization make this inappropriate (e. g. using different scaling factors for each gene or anything like that)?

I have tried to read the TMM-normalization paper (cited below) and it seems that when there are more than two samples, it picks a reference sample and calculates a single scaling factor for each pairwise comparison between the reference and each non-reference sample. That would seem to suggest that the within-sample comparison like the one outlined above is possible. Is it?

Robinson, M. and A. Oshlack (2010). "A scaling normalization method for differential expression analysis of RNA-seq data." Genome Biology 11(3): R25.

Answer: TMM won't change the order of values, it's meant to adjust between samples only. So if you're asking, "Does gene A have a higher FPKM than gene B in sample 3", then you'll get the same answer with/without normalization (the exact size of the difference might differ, though).
Question: I have a dataframe containing a list of enriched GO terms resulting from a differential expression analysis. I want to filter out GO terms of interest by filtering with a second dataframe containing those GO terms I'm looking for.

Based on this post

[https://community.rstudio.com/t/dplyr-filter-from-another-dataframe/7207][1]

I wrote the following code

Extracted_GO <- filter(BP_1_5, GO_Accession_Number == GO_Search$GO_Accession_Number)

which gives the following error

longer object length is not a multiple of shorter object length

and

str(BP_1_5)

gives

tibble [5,325 x 5] (S3: tbl_df/tbl/data.frame)
$ GO_Accession_Number: chr [1:5325] "GO:0000028" "GO:0000045" "GO:0000045" "GO:0000070" ...
$ Biological.Process : chr [1:5325] "ribosomal small subunit assembly" "autophagosome assembly" "autophagosome assembly" "mitotic sister chromatid segregation" ...
$ HGNC_Symbol : chr [1:5325] "ERAL1" "UBQLN1" "STX12" "INCENP" ...
$ Gene.name : chr [1:5325] "Eral1" "Ubqln1" "Stx12" "Incenp" ...
$ GO.domain : chr [1:5325] "P" "P" "P" "P" ... data.frame':	5325 obs. of 6 variables:

and

str(GO_Search)

gives

tibble [1,727 x 2] (S3: tbl_df/tbl/data.frame)
$ GO_Accession_Number: chr [1:1727] "GO:0062201" "GO:0098610" "GO:0075003" "GO:0075002" ...
$ Description : chr [1:1727] "actin wave" "adhesion between unicellular organisms" "adhesion of symbiont appressorium to host" "adhesion of symbiont germination tube to host" ...

Given my two dataframes are of fixed length and can't be changed is there a way of achieving my goal?

[1]: https://community.rstudio.com/t/dplyr-filter-from-another-dataframe/7207

Answer: filtered <- BP_1_5[BP_1_5$GO_Accession_Number %in% GO_Search$GO_Accession_Number, ]
Question: Hello,

I would like to ask if you know any R package which will help me in doing a mass search for genes relations and which pathways are they corresponding to?
E.g I had about 34 thousand genes (i.e ENSG00000146648 - coded in ensembl DB) derived from GBM cancer patients. I want to list their functions and check in which pathways each of them is involved. Do you have any ideas?

I was looking for a similar question, but I could not find one (expect biopython, which I am not using - unless I will be forced to. ;))

Answer: You can get GO annotations for genes in Ensembl. Use either the [Ensembl perl API][1] or the Bioconductor [biomaRt R package][2].

[1]: https://www.ensembl.org/info/docs/api/index.html
[2]: http://bioconductor.org/packages/release/bioc/html/biomaRt.html
Question: Hello,

After running Varscan and Mutect on a set of 10 patients (tumor / normal comparison), I have run through a pipeline of false-positive filtering. When I look at my resulting Ts/Tv ratio (by manual calculation, snpEff summary file, SnpSift tstv calculation or GATK VariantEval), it is quite low for human whole genome sequence data (1.3-1.6). I have read all I can find here and in papers about the expected ratio, and how a low ratio could denote a great deal of false positives.

I ran Varscan with relatively lax parameters for calling somatic mutations (5 reads in N, 8 in T, but strand bias filtered), however I thought Mutect would call a confident set. Both SNP callers end up with a low Ts/Tv. My question is, can I chalk this result up to false positives (which is okay with me, I wanted a sensitive not specific call set), or could it be a problem with the BAM alignment? I suppose a poorly aligned BAM would lead to false positives too, but any insight or information would be greatly appreciated.

Answer: We should expect Ts/Tv ratio of somatic point mutations to be wildly variable across tumor types... depending on various mutagens, or the mechanisms involved in DNA repair. I can't seem to find a publication that confirms this assumption, but [this figure comes close][1]. Here are my quick and dirty Ts/Tv ratios of mutation calls grabbed from [that paper][2], but please double-check my work.

[1]: http://www.nature.com/nature/journal/v500/n7463/fig_tab/nature12477_F2.html
[2]: http://www.nature.com/nature/journal/v500/n7463/full/nature12477.html

**Note**: A caveat in the data below is that some cohorts are exomes while others are whole-genomes. Since there's more GC content in exomes, these Ts/Tv ratios are not perfectly comparable... but good enough for our point to hold.

```
Cancer Type Ts/Tv
ALL 0.949906
AML 2.128909
Bladder 1.325778
Breast 0.859808
Cervix 1.265049
CLL 1.006487
Colorectum 2.163191
Esophageal 1.38155
Glioblastoma 3.53876
Glioma Low Grade 2.244252
Head and Neck 1.172555
Kidney Chromophobe 2.545455
Kidney Clear Cell 1.165541
Kidney Papillary 1.116037
Liver 1.222369
Lung Adeno 0.439277
Lung Small Cell 0.569885
Lung Squamous 0.635106
Lymphoma B-cell 0.971431
Medulloblastoma 1.381825
Melanoma 8.54497
Myeloma 1.303654
Neuroblastoma 0.566366
Ovary 0.876746
Pancreas 1.021448
Pilocytic Astrocytoma 1.837178
Prostate 1.220668
Stomach 3.006267
Thyroid 2.161623
Uterus 1.632635
```
Question: Hello,
I am doing a RNA-Seq analysis with R using the limma package. I have RNA-Seq data of different mutant lines of a model organism which have two different phenotypes in comparison with the siblings as control. The samples are like that:

- line 1, phenotype 1, mutant and sibling
- line 2, phenotype 1, mutant and sibling
- line 3, phenotype 2, mutant and sibling
- line 4, phenotype 2, mutant and sibling
- line 5, phenotype 2, mutant and sibling

That means line 1 and 2 have the similar phenotype 1
and line 3, 4, and 5 have the similar phenotype 2
which differ from phenotype 1. Let's say phenotype 1 has more cells than the control and phenotype 2 has less cells than the control.

I did the RNA-Seq analysis for each line, but now I want to compare the lines of one phenotype. I want to find the genes which are different in the lines with the same phenotype in comparison to the siblings. My makeContrast command looks like that:

cont.matrix <- makeContrasts(pheno1 = (line1.mut+line2.mut)-(line1.sib+line2.sib),
pheno2 = (line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib),
levels=design)

Finally I want to find genes (eg. genes involved in the cell cycle) which are in phenotype 1 up and in phenotype 2 down regulated or the other way around.

cont.matrix <- makeContrasts(pheno1vspheno2 = ((line1.mut+line2.mut)-(line1.sib+line2.sib)) -
((line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib)),
levels=design)

I am not sure if a can do a contrast matrix like this. Do you think this will give me the genes I am interested in?

Answer: You may want to double-check this but it seems to me that with:

pheno1vspheno2 = ((line1.mut+line2.mut)-(line1.sib+line2.sib)) -
((line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib))

You are looking for an interaction between mutant/sibling state and phenotype. *I.e.* genes that respond differently between mutant and sibling depending on whether they are in phenotype 1 or 2.

**EDIT** russhh in his nice answer points out that the contrasts should be averaged. So the above should in fact be (please check ok):

pheno1vspheno2 = ((line1.mut+line2.mut)/2 - (line1.sib+line2.sib)/2) -
((line3.mut+line4.mut+line5.mut)/3 - (line3.sib+line4.sib+line5.sib)/3)

----

If you want the difference between phenotype 1 and 2 regardless of mutant/sibling state you could use:

pheno1vspheno2 = (line1.mut + line2.mut + line1.sib + line2.sib)/4 -
(line3.mut + line4.mut + line5.mut + line3.sib + line4.sib + line5.sib)/6,

I would suggest you add to your count matrix one or more a dummy genes with behaviour that you want to pick up as differential or not significant to test whether the contrasts do what you expect.
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: Flybase has direct access to its postgres chado database.

http://flybase.org/forums/viewtopic.php?f=14&t=114

hostname: flybase.org
port: 5432
username: flybase
password: no password
database name: flybase

e.g.

psql -h flybase.org -U flybase flybase
Question: Good afternoon,

I have a question about using collapseReplicates in DESeq2.
As far as I understand, this function adds up the counts belonging to one biosample. I understand the meaning of this if the number of technical replicates per biosample is the same for all samples. Please tell me what should I do if the number of technical replicates per sample differs?

For example, SAMNXXXXX corresponds to SRRXXXXX1 (count = 5) and SRRXXXXX2 (count = 6), SAMNYYYYY corresponds only to SRRYYYYY1 (count = 10). If I add up the counts for SAMNXXXXX (5 + 6 = 11) and then compare it with count for SAMNYYYYY (10), I will get an incorrect conclusion that the expression is higher in SAMNXXXXX.

Maybe I need to take the arithmetic mean or something else? It seems to me that the arithmetic mean is not very reasonable. For example, I have counts of 181 and 2 for different replicates of the same biosample.

Note: this situation is not observed for most samples. For example, in a particular dataset there are 89 biosamples without technical replicates and and 5 biosamples with 2 technical replicates in each.

Thanks!

Good regards,
Poecile

Answer: Collasping **technical** replicates takes place before DESeq2 internal normalization, so you don't need to worry about the arithmetrics. Difference in sequencing depth will be handled exactly as with non-collapsed samples: by dividing the counts by a size factor calculated based on the median of ratio method. Therefore, it is ok to not have the same number of technical replicates by biosample.
Question: Hello everyone

I am doing network analyses by using WGCNA, I got error message when I am setting the soft threshold. Here is what I got.

I am following these commands:

> library(WGCNA)

disableWGCNAThreads()
> library(cluster)
> options(stringsAsFactors = FALSE)
> femData = read.csv("gene.csv")
> names(femData)
[1] "gene_id" "gene" "locus" "log2.fold_change."
[5] "test_stat" "p_value" "q_value" "AF"
[9] "AF.1" "AF.2" "AD" "AD.1"
[13] "AD.2"
> dim(femData)
[1] 6765 13
> datExprFemale=as.data.frame(t(femData[, -c(1:7)]))
> names(datExprFemale)=femData$gene
> rownames(datExprFemale)=names(femData)[-c(1:7)]
> powers=c(1:10)
> sft=pickSoftThreshold(datExprFemale,powerVector=powers)
Error in summary(lm1)$coefficients[2, 1] : subscript out of bounds

In addition: Warning message:

executing %dopar% sequentially: no parallel backend registered

I tried to find similar questions in archive, but I could not.

Please I need help with this issue if anyone able to help.

Best

Shaima

Answer: The error means you're trying to access an element outside the array; this usually means one your indices is greater than the last available index of the array.

As for the warning, check the [doParallel docs][1]: "Remember: unless registerDoMC is called, foreach will not run in parallel. Simply loading the doParallel package is not enough."

[1]: https://cran.r-project.org/web/packages/doParallel/vignettes/gettingstartedParallel.pdf
Question: Hi,

I would appreciate it if anyone could direct me on how to give different colors for the up and down-regulated genes using the R Bioconductor package **enhancedvolcano**.

Just as indicated in the following link.

Thanks!

https://galaxyproject.github.io/training-material/topics/transcriptomics/tutorials/rna-seq-viz-with-volcanoplot/tutorial.html#:~:text=A%20volcano%20plot%20is%20a,the%20most%20biologically%20significant%20genes.

Answer: Hey,

This can be achieved via customised coding or EnhancedVolcano. Please take a look at the vignette, where there is a dedicated section for this (see Section 4.9): http://bioconductor.org/packages/devel/bioc/vignettes/EnhancedVolcano/inst/doc/EnhancedVolcano.html

Kevin
Question: Hi,

I have some paired end exome sequencing, and after checking out the samples in fastqc, I found there was an enrichment for Illumina adapters. I used Trimmomatic to remove the adapters (using Illumina Clip), which gives me paired reads, and unpaired reads for forward and reverse respectively post trimming.

I'm following the GATK best practises for variant calling, which recommends BWA MEM for alignment. I can use the paired, trimmed reads for to align, but there's quite a large proportion of reads which are trimmed and unpaired. The bedrock of my question is how can I use these effectively?

Would the best way be to do an alignment with the pairs, then unpaired, and combine the SAM files in some way? Or even merge the BAMs?

Anyone had a similar experience? Advice is welcome!

Thanks

Answer: I do this regularly with exome datasets.

Reads where one mate fails but are otherwise fine are included in a single-end file. It's still good data! Map them the same way you do the PE data, and subsequently combine the BAMs using Picard's MergeSamFiles or another approach.

Some cleaning approaches make this easy. I use expHTS to manage my cleaning now - does a fantastic job and is blazingly fast with stream-based processing. Documentation is forthcoming, but the developers are happy to help.
Question: Dear all,

I was recently asked in a lab meeting whether there may be a gene length bias in the results produced by WGCNA, if the input data consists of normalized counts or variance-stabilized transformed data from DESeq2 (which does not correct for gene length). I am not sure how to answer this, as [there are several papers that used DESeq2 normalized counts as input for WGCNA][1], and I thought that this was actually [recommended by the authors of WGCNA][2].

In other words, could my detected expression modules be biased for gene length, with some modules being particularly driven by the length of the genes in that module as opposed to its actual expression?

[1]: https://scholar.google.com/scholar?hl=en&as_sdt=0%2C9&q=deseq2%20wgcna&btnG=
[2]: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html

Answer: In any analysis that is based on the number of reads mapping to a gene, the longer a gene is the more reads will map to it.

In DE analysis (such as DESeq and edgeR), this means that genes that a long gene is more likely to be called significantly differential than a short gene as there is less noise.

I don't know of any literature on how WGCNA might be affected by this, however, one can imagine that where the expression is more accurately estimated (as it is with longer genes) there is a higher chance of the correlation being more significant. So, yes, I might image that you would find that longer genes are more likely to be assigned to clusters.

Might even be a short paper in it if you could demonstrate it.
Question: Hi, I was running plink(1.9) tool to create bed file (`plink --file inputfile --make-bed --noweb --out outputfile`) which worked fine, but when I switched to plink2, the same command doesn't work. Can we do everything that plink(1.9) does with plink2? Can someone please clarify this? Thank you.

Answer: Sorry about not noticing this earlier.

The main difference is that plink 1.9 is essentially finished, while plink 2.0 is an alpha-stage program which will have significant unfinished components for a while to come. As a consequence, current development priorities for plink 2.0 are centered around things which are impossible to do with plink 1.9, such as handling multiallelic/phased variants and dosage data and reliably tracking REF/ALT alleles; while things that plink 1.9 already handles perfectly well, such as working with .ped/.map file pairs, have been deliberately deprioritized for now.

So, you should stick to 1.9 as long as it's good enough for the jobs you need to perform. But once you need to do something outside 1.9's scope, you're likely to find that 2.0 already has the additional feature you need (or it'll be added quickly after you ask for it).