cannin/biostars_qa huggingface dataset:
Question: Does anyone know if it is possible to specify which assembly to use when constructing a query for Entrez ?

For example, if I do such a query with EDirect:

esearch -db gene -query "brca1 [ALL]human[ORGN]" -sort "relevance" | \
efetch -format docsum | \
xtract -pattern DocumentSummary -element Name MapLocation Description OtherAliases Id \
-block GenomicInfo -element ChrLoc ChrAccVer ChrStart ChrStop

the GenomicInfo that I get seems to be according to the GRCh38 assembly. The application that I'm developing needs to use the GRCh37 assembly, however.

Any help is much appreciated.

Answer: Actually, it turns out that this can be done. I wrote to NCBI and got a response.

The key is to look in the **LocationHistType** block for:

- a specific annotation release (see here and here for some explanation of annotation releases). For example, GRCh37.p13 is coded by NCBI as Annotation Release *105*.

- the corresponding assembly accession, which is a RefSeq Assembly ID. For GRCh37.p13 it is *GCF_000001405.25* (see here for more information).

The EDirect commands should look something like this:

esearch -db gene -query "brca1 [ALL] AND human [ORGN]" | \
efetch -format docsum | \
xtract -pattern DocumentSummary \
-element Name MapLocation Description OtherAliases Id ChrLoc \
-block LocationHistType \
-match "AnnotationRelease:105" -and "AssemblyAccVer:GCF_000001405.25" \
-element ChrAccVer ChrStart ChrStop
Question: Hello to all and I hope you are all doing well. I am a biology student who's doing a study on bird gut microbiota analysis thorugh 16S rRNA. I collected 3 biological samples per time point (we have five different time points) and I want to see the changes in bird's gut microbiota across those five time points. However, one of the time points only has two biological replicates as the other birds did not have gut contents for me to extract. I have extra birds per time point in case anything happens on the original biological replicates. However, for this particular time point, all the extra bird's and one of the replicate's gut were empty so I'm stuck with only two workable biological replicates for this time point. Now I am having a concern since most microbiome studies use a minimum of biological triplicates and I've read in many statistical reviews that higher biological replicates is necessary to reliably infer the average or variation in a microbial population. I'm wondering if there's a way for me to statistically justify this discrepancy in biological replicates or a way for me to proceed the analysis of this one, or is there any feasible way for me to analyze samples with uneven biological replicates? Thank you so much and apologies if I'm using statistical terms wrong as I am not completely versed in the field. Thanks again and regards.

Answer: The good news is that this is a common problem, and if you have no way of obtaining additional data, I don't think you'd be faulted for just analyzing what you have. Review may be difficult if you're looking to publish formally, but you should still analyze it and explain what happened.

Simple methods like t-tests can still be applied for comparison between n=3 and n=2 groups. However, there are more correct ways of handling microbiome missing data and timeseries analysis than doing pairwise comparisons. The bad news is that these methods are generally high-level and require some programming knowledge to apply. An example would be [this publication][1], in which they discuss the problem you're having with missing data. The software they use can be found [here][2].

In general, microbiome data (and 16S in particular) are very noisy, and it often takes very large samples sizes to make biologically-relevant inferences. It is great that you have done timeseries though, as it will improve your statistical power by a lot. Typically, there is high variability between individuals' microbiomes, but collecting data over time lets you model the variability within each individual's microbiome as well as between them. Probably the best example of a high-powered microbiome time-series study is the [HMP2 publication][3], but even with a very large sample size, they weren't able to say much biologically speaking. So don't be disappointed if the analysis doesn't find much (my first 16S study didn't); it is still valuable to analyze and share it, since it contributes to our global database of 16S data.

[1]: https://academic.oup.com/bioinformatics/article/34/3/372/4157442
[2]: https://github.com/tare/GPMicrobiome
[3]: https://www.nature.com/articles/s41586-019-1237-9
Question: Hi Biostars:

I am facing a problem with differential expression analysis, were due to the intrinsic features of the samples we can't have replicates in one condition.

**Study design:** We have a patient with a very rare translocation , no other similar translocation has been described in the world and we expect that expression of the genes surrounding the translocation is altered. Hence, we have perform RNA-seq of the translocated patient and 4 controls.

**How should I proceed?**:

Before you kill me for asking "Can I do differential expression without replicates?", I known that [EdgeR][1] and [DEseq2][2] provide ways to proceed without replicates , and that [NOISeq][3] could be used without replicates. However, in this case we have replicates for the controls were the biological variance of each gene could be estimated, but no for the the patient, as there is no other individual in the world. So, **my question is:** Is there any way of estimating the variance for the control group, and then compare to the expression in one single sample, the patient in this case? Or better, do a dispersion estimation for the translocation patient?

**my options** (From EdgeR docs):

- Use the genes and transcripts that are far away or in other chromosomes than the chromosome with the translocation to estimate the dispersion of that sample
- Use a dispersion value defined previously.

Have any of you faced a similar problem, and, furthermore, have anybody tested how do the two options above mentioned that EdgeR provide for working without replicates perform?

thanks for reading :)

[1]: https://www.bioconductor.org/packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
[2]: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
[3]: https://www.bioconductor.org/packages/release/bioc/vignettes/NOISeq/inst/doc/NOISeq.pdf

Answer: The short answer is that you just proceed as usual. limma, edgeR and DESeq2 have no trouble with this scenario, although the edgeR quasi-likelihood pipeline would be better than the other options. The packages simply estimate variability from groups where you do have replication (controls in your case), and apply the same dispersion estimates to all the samples in all the groups.

The edgeR and DESeq2 pipelines for no replicates are for when *none* of the groups have any replicates. You however do have replicate controls.

edgeR can be used right down to a two-group comparison with n=2 in one group and n=1 in the other. I'm not saying that such small sample sizes are desirable, but the package will do the best it can with what it gets and will present scientifically defensible results even in that extreme scenario. You can see an example of an n=2 vs n=1 analysis in the discussion to this paper (i.e., my reply to Conrad Burden's first report):
https://f1000research.com/articles/5-1438

BTW, the same question has been asked several times on the Bioconductor Support forum, for example: https://support.bioconductor.org/p/63585/ or https://support.bioconductor.org/p/61904/
Question: Hi :)

A few weeks ago I was looking for a tool that would help me get "DNA composition" statistics for my sequencing. Something that would give me a dataset which which I could ask questions about GC bias, or over-represented sequences, motifs, etc. There are tools to answer each question specifically, but I was looking for something more general from which many analyses could be built on. This led me to k-mer counting, and all the down-stream tools which leverage k-mer result files.

k-mer counting tools are pretty cool, but all the ones I tried had some draw backs; like high RAM requirements, very long run-times (although the latests bloom-filter based tools seem to mitigate this somewhat), but most importantly requiring a specific k-mer size to be chosen. I really wanted all mers in the dataset so I could look at ''GC' and 'GCG' and 'GCCGACGGACGAC' without having to re-run any analyses. I couldn't find a tool like this after a brief search, so gave up and wrote my own in NumPy based off suffix-arrays.

Two weeks later I have a functional program in the sense that I get results, but before I invest any time making it usable for others, I thought I should investigate further if there are tools which already do this. Making a suffix array was a nice learning experience for me so I haven't lost anything if such tools already exist - and if they do i'd love to compare performance characteristics - but if not I might consider tidying up the code and making proper documentation. Does anyone know of such tools?

Thank you so much, and happy Diwali :)

Answer: This is an interesting question because even though there are dozens of k-mer tools available they do all seem to have the same limitations: very short k-mer length requirements (usually <32bp) and fixed length search. It seems like there are some new tools reporting performance improvements but they still seem to have those features. Tallymer (part of GenomeTools) is the only tool, to my knowledge, that can compute statistics for an arbitrary length (or range) of k-mers. It sounds like you want the 'occratio' subcommand, which runs for a range of k-mers (as opposed to the 'mkindex' command which runs for a fixed length). Here is an example to get a range of unique/nonunique k-mers and get the relative amounts.

gt tallymer occratio -output unique nonunique -minmersize 10 -maxmersize 180 -esa db
gt tallymer occratio -output unique relative -minmersize 10 -maxmersize 180 -esa db

Note, you have to create 'db' first, but this was just an example. I can't go into the details, but I know that Tallymer uses a data structure called an enhanced suffix array. This is by the same author of Vmatch, and MUMmer (and many other tools). I should mention that Vmatch is also amazingly flexible for doing k-mer arithmetic and we could probably come up with an equivalent command using that tool. I agree with you about learning by doing, but this is a problem I would personally not try to solve since there are great tools available that allow you to answer virtually any question you may have about genomes and k-mers.
Question: Have anyone use Araport11 genome annotation file to replace TAIR10 genome annotation file?Because Araport11's GTF was released on 2016 ,which is newer than TAIR's GTF. And Araport11 also provide gene name,like "gene-id:AT1G01060; gene name"LHY" " I want to ask whether Araport11 is better for use as genome annotation file?

Answer: that would be up to you to decide. The annotation provide by Araport11 is indeed newer than Tair10. They should be nonetheless interchangeable, however general rule is to use the one from which you also use the genomic sequence (though the genome for TAIR10 and Araport11 should be the same).

If Araport 11 suits your needs better then you can safely go for that one.
Question: samtools flagstat results:

some number + 0 in total (QC-passed reads + QC-failed reads)

my this "some number" does not match actual number of reads in my paired fastq file. Is this expected?

Answer: The number of reads reported is the number actually in the file. If your aligner produced secondary alignments then this will often be higher than the original number in the fastq files.
Question: Here is an example of volcano plot produced using edgeR results on DEG analysis. So my question is why there is no values in the middle of plot (indicated by red circle)?

![Volcano plot][1]

[1]: https://i.pinimg.com/originals/07/69/3d/07693dfcb3c8fe55fcae39d7caa64239.jpg

Answer: This is because there is a relationship between p-value and fold change. With no fold change (log2FC = 0) you have no significant difference hence p-value of 1 which is zero after -10log. With higher FC you can get lower p-values (hence higher -log10 values).
Question: Hi, there

I downloaded RSEM expected_count data of TCGA TARGET GTEX from UCSC Xena to perform differential expressed analysis. I noticed the unit is log2(x+1), but I failed to turn the data into integers. The input of DESeq2 should be a matrix with integers. Any suggestions? Thanks in advance!

Best,

Issac

Answer: log2(x+1) data is not suitable for analysis with DESeq2 (nor edgeR/voom). So you would at least have to convert it back to counts by doing (2^x) -1. You could then simply round to the nearest integer. Alternatively you can download un-tarnsformed RSEM output from Firehose at https://gdac.broadinstitute.org/
Question: Dear all,

I am new to the field. I am trying to analyze single end 100b FastQ files with ~70million reads/sample. I am trying to determine if adapter sequences are present and if so how to go about them. I ran FastQC on the files and reports show they each have an "overrepresented sequence" of an "illumina index adapter" in them.

I have the following questions:

1. Does sample1 look like a trimmed file or it requires adapter trimming?

2. If further trimming is recommended what would be the best seq/adapter option to be used for cutadapt/TrimGalore? [See below for my thoughts so far]

3. Based on the FastQC report, do I need to worry about presence of any other adapter sequences beside the index?

My thoughts on question 2:
The sequences for illumina index adapter format appear to be:

GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG

These are the adapter sequences found in my FastQC report for sample 1:

GATCGGAAGAGCACACGTCTGAACTCCAGTCACCATGGCATCTCGTATGC
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACCATGGCATCTCGTATG

I am thinking of using below options for cutadapt/trimgalore to remove the adapter(s):

trim_galore sample1.fastq.gz -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -q 20 --length 20 –fastqc

However, it seems that trimmomatics for instance only takes care of the initial sequence of the index adapter (only up to Ns and not after):
https://github.com/timflutre/trimmomatic/blob/master/adapters/TruSeq3-SE.fa

Many thanks for your time and reply beforehand.

Answer: The sequence to trim if having Universal Adaprer contaminations is `AGATCGGAAGAGC`.
In you case (0.5%) I would not even bother and directly align the files without any manipulations.
Question: Hi,

I'm new to discoSNP (using 2.2.1).

It seems that the script simple_test.sh was not designed for this version.

the command line using fof.txt works fine.

I think you should remove simple_test.sh from the package., it's confusing.

Answer: Hi Sebastien,

Thanks for this message. Indeed, simple_test should not be part of the release. It won't appear in the 2.2.2 release.

Pierre
Question: Hellow I have a question on how to calculate fold changes when analyzing gene expression changes between multiple tumor and control samples per gene?

Answer: Or the bioconductor limma package if you are dealing with arrays and/or RNA-Seq to analyze your data

Limma will give you the log2 expression changes based upon statistical values
Question: I have five SRR files for one sample (five runs, Single read - Illumina experiment). What are the steps to be taken before further analysis? I am converting this files to fastq format. Should I merge this five fastq files before mapping to the genome or map them separately. I am new to NGS and don't fully understand what "runs" means. Should I understand them as five independent experiment repeats and analyse separately?

Answer: Hi
Most mappers accept multiple fastq files, so no need of merging them as long as read1 and read2 files are properly supplied. It is also best to retain the original header line of the fastq (-F argument to fastq-dump). See here https://edwards.sdsu.edu/research/fastq-dump/ . The header can be important for FastQC and other downstream applications. SRA classifies as Sample (SRS/ERS) -> Experiment (SRX/ERX) ->Data (SRR/ERR). SRR* files are the resulting data files of a particular experiment (SRX) run on a particular sample (SRP). Running same experiment on the same sample multiple times can result in multiple SRR files.
Question: Hi, i want to phase my data in vcf format using impute2's "-phase" option. So i tried to convert vcf file with plink:

plink --vcf data.vcf --recode oxford --out converted_data

But produced .map file does not conform to the required three-column format.

Is there any other way to convert vcf to impute2 format?

Answer: Can you clarify what the problem is? `--recode oxford` should produce a .gen file and a .sample file; this command shouldn't generate a .map file at all.
Question: I have RNAseq data on B. subtilis, and mapped it to the genome so I have gene names and even BSU numbers for all of the genes. So far so good. Now I am struggling to find and assign functional annotations to each, preferably GO annotation.

Despite B. subtilis being a model organism, it is not included in the GO consortium database; they only have E. coli.

KEGG and bsubcyc both have B. subtilis functional annotations, but I cannot find a way to download any flat file from either site that associates each gene with its functional assignments.

How do people do this? What set of bioinformatics skill am I totally unaware of but that seem to be taken for granted? Do people "webpage scrape" for this info from KEGG? Is there a database format I need to learn?

Answer: The most straight-forward approach to obtaining GO annotation would in my opinion be to fetch it from UniProt. You can do that in several ways, for example, by querying for it via their REST API:

http://www.uniprot.org/uniprot/?query=taxonomy:224308&columns=id,go-id&format=tab
Question: I made 3 binary files from MAP and PED files.

I'd like to know how to get individual genotypes per specific SNP.

When I open the .bim file, I just see the SNP and allele information but individual information (i.e. sample Id)

How to map the information?

I'd like to get information like this, and could open in Python or R.

Do I need some other tools for that? please let me know. For example,

SNP SAMPLE GENOTYPE
rs53576 1111 AA
rs53576 1112 GA
rs53576 1113 GG
rs53576 1114 GA

Thank you in advance.

Answer: plink --bfile --extract --recode --out

This should generate the desired output
Question: Hi everybody,

can anyone guide me from where to download the hs37d5.fa genome merged with phiX genome?

I have downloaded few BAM files containing the tag in header section @SQ SN:PhiX LN:6386 together with other hs37d5/fa contigs. I am looking for which version of the reference sequence was used for these files.

any help would be appreciated!

Regards

(Note: I am aware of phiX genome is used perhaps a control in illumina sequencing.)

Answer: If you want *exactly* the MPI-EVA version (which includes a circularised phiX and extended reference mtDNA) the construction is similar to the reference used by the 1000 Genomes Project (compare ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/),
with minor changes. It's made as follows :

1 Download individual chrs from ensembl ftp (just like 1000g)
ftp://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/

2a Download the newer version of the mitochondrion (NC_012920, just like 1000g)
http://www.ncbi.nlm.nih.gov/nuccore/251831106

2b Copy the first 1000bp of the mitochondrion onto its end. The resulting sequence is named "MT".

3 Download the concatenated decoy sequences from 1000 Genomes:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5cs.fa.gz

Also compare their READMEs:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/README_human_reference_20110707
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.slides.pdf

4 Download the Human herpes virus (NC_007605, aka EBV) from NCBI, just like 1000g. The sequence is then named "NC_007605".
http://www.ncbi.nlm.nih.gov/nuccore/NC_007605

5a Download phiX-174 reference (NC_001422).
http://www.ncbi.nlm.nih.gov/nuccore/NC_001422

5b Copy the first 1000bp of phiX onto its end, name the result "phiX".

6 Create a reference (whole_genome.fa) with chrs 1-22, X, Y, extended
NC_012920 MT, the non-chromosomal supercontigs, the
NC_007605 EBV, the decoy sequences (hs37d5), extended phiX.
The order is chosen to match 1000 Genomes (plus phiX), see their fai file:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz.fai

Note that two sequences (MT, PhiX) are circular and have been extended to
facilate alignment. The correct incantation to wrap these alignments to their
correct length is

bam-rewrap MT:16569 phiX:5386

or

bam-rmdup -z MT:16569 -z phiX:5386
Question: Hi there,

When I am running the DESeq pipeline on the dds object and getting this error message.

> dds_res<-DESeq(dds_PvsN)
estimating size factors
Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc, :
every gene contains at least one zero, cannot compute log geometric means
In addition: Warning message:
In class(object) <- "environment" :
Setting class(x) to "environment" sets attribute to NULL; result will no longer be an S4 object

What strategy should be applied to resolve this conflict?

Thank you very much!

Imran

Answer: I encountered this error only once in the past. It is as stated: every gene in your data has at least one zero value, and this creates an issue for the size-factor calculation.

Solutions:

1. add a pseudo-count value of '1' to your data
2. use: `estimateSizeFactors(dds_PvsN, type = 'iterate')`

Kevin
Question: I want to gseGO() provided by clusterProfiler R package. My gene list is the unregulated markers from a cluster generated by UMAP. The number is log2 fold change.

```r
str(head(cluster7.genes))
Named num [1:6] 2.14 1.81 1.76 1.71 1.67 ...
- attr(*, "names")= chr [1:6] "10082" "9353" "53353" "55553" ...

cluster7.enriched <- gseGO(geneList = cluster7.genes, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pvalueCutoff = 1, verbose = FALSE)

Warning messages:
1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
There are duplicate gene names, fgsea may produce unexpected results.
2: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
All values in the stats vector are greater than zero and scoreType is "std", maybe you should switch to scoreType = "pos".
```

But I don't know why it returns errors even when the `pvalueCutoff` is 1. Do I need to change `scoreType` to `"pos"` to solve this error? If so, how to do it since I can't see any relevant parameters available in `gseGO()` for me to change? Thanks.

![enter image description here][1]

[1]: /media/images/b5d2edcb-d3a1-4e12-b5e3-ed40e6e3

Answer: It has given you a warning message, not an error. You should take a look at the contents of cluster7.enriched to see if it has yielded any results.

When doing a preranked GSEA like this, the inputs are generally a mixture of up and down regulated genes, so there will be positive and negative numbers. Some folks like to use the log fold changes as input. I like to use the DESeq test statistic.

Also it looks as though you have duplicate gene names. This could be a problem and it is best to explicitly resolve this before GSEA. Some approaches used are to take the mean of duplicate genes, or just keep the entry with the most extreme value.
Question: Hi all.

Does any know of a script or code that can convert a GVF file to VCF? There is a reference to a link https://code.google.com/p/gvf2vcf/ but it is not accessible anymore.

Any help will be appreciated.

Thanks

Answer: Found a perl script here that should do it (scroll to bottom for usage instructions)