cannin/biostars_qa huggingface dataset:
Question: Hi,

I'm using mirdeep2 for finding miRNAs on mouse data.

I ran mapper.pl as a first step. It runs well and give me output. But throws an error at the end about mapping statistics. Does anyone know what this error is? And I'm not sure if I should rely on my results because of this error.

I ran using the command:

mapper.pl Sample.fastq -e -h -p ~/refs/mm9 -s Sample.fa -t Sample.arf

And I get this error at the end:

# reads processed: 24468373
# reads with at least one reported alignment: 13415055 (54.83%)
# reads that failed to align: 10102035 (41.29%)
# reads with alignments suppressed due to -m: 951283 (3.89%)
Reported 18449906 alignments to 1 output stream(s)
Mapping statistics

#desc total mapped unmapped %mapped %unmapped
Use of uninitialized value $count2 in subtraction (-) at /opt/sam/miRDeep/2.0.5/mapper.pl line 630, line 18449906.
Use of uninitialized value $count in subtraction (-) at /opt/sam/miRDeep/2.0.5/mapper.pl line 630, line 18449906.
total: Use of uninitialized value $count in print at /opt/sam/miRDeep/2.0.5/mapper.pl line 630, line 18449906.
Use of uninitialized value $count2 in print at /opt/sam/miRDeep/2.0.5/mapper.pl line 630, line 18449906.
0 Use of uninitialized value $count in division (/) at /opt/sam/miRDeep/2.0.5/mapper.pl line 631, line 18449906.
Use of uninitialized value $count2 in division (/) at /opt/sam/miRDeep/2.0.5/mapper.pl line 631, line 18449906.
Illegal division by zero at /opt/sam/miRDeep/2.0.5/mapper.pl line 631, line 18449906.

Answer: I was able to solve the issue. I included a new option `-m`, and it started working well.

mapper.pl sample.fastq -e -h -m -v -p ~/refs/mm9 -s Sample.fa -t Sample.arf

Results:

parsing fastq to fasta format
discarding short reads
collapsing reads
mapping reads to genome index
# reads processed: 6409517
# reads with at least one reported alignment: 4426074 (69.05%)
# reads that failed to align: 1850959 (28.88%)
# reads with alignments suppressed due to -m: 132484 (2.07%)

Reported 4984883 alignments to 1 output stream(s)
trimming unmapped nts in the 3' ends
Mapping statistics
#desc total mapped unmapped %mapped %unmapped
total: 17205138 9753961 7451177 0.567 0.433
seq: 17205138 9753961 7451177 0.567 0.433
Question: I have a list of protein GI and would like to get accession number of the genome (DBSOURCE in genbank file) using Eutilities.

What would be the easiest way to get genome accession numbers for a list of protein GI?

Answer: For future reference:

There is an interesting file @NCBI which contains information about various accession numbers/gene names: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2accession.gz

This file is automatically generated everyday and should contain latest information about GeneID's and accession numbers.
Question: I have a GWAS sumstats file called `sumstats.txt` in the following format:

SNP Freq.A1 CHR BP A1 A2 OR SE P
10:100968448:T:AA 0.3519 10 100968448 t aa 1.0024 0.01 0.812
10:101574552:A:ATG 0.4493 10 101574552 a atg 0.98906 0.0097 0.2585
10:10222597:AT:A 0 10 10222597 a at 0.9997 0.01 0.9777

For every entry in the sumstats file I need to munge column 1 to retain only the chromosomal information, effectively deleting everything after the first colon, whilst leaving the rest of the file intact.

The final file would look like this:

SNP Freq.A1 CHR BP A1 A2 OR SE P
10 0.3519 10 100968448 t aa 1.0024 0.01 0.812
10 0.4493 10 101574552 a atg 0.98906 0.0097 0.2585
10 0 10 10222597 a at 0.9997 0.01 0.9777

My `awk`/`sed` knowledge is not great but I've had a bash (npi!!):

echo "10:100968448:T:AA 0.3519 10 100968448 t aa 1.0024 0.01 0.812" | awk '{n=split($0,a,/[:_]/); print a[1]"\t"a[2]"\t"a[2]+1"\t"a[3]"/"a[4];}'

Giving:

10 100968448 100968449 T/AA 0.3519 10 100968448 t aa 1.0024 0.01 0.812

But not sure how to get rid of the info after the colon and pipe the file in correctly. Also tried:

echo "10:100968448:T:AA 0.3519 10 100968448 t aa 1.0024 0.01 0.812" | awk '{n=split($0,a,/[:_]/); print a[1];}'

Returning:

10

I'm looking for a `bash`/`awk`/`sed` one liner for this as I have quite a few sumstats files to process.

Any suggestions would be greatly appreciated.

Answer: It looks like your goal is to split the first field on the colon, then print the rest of the file. Your awk script is splitting the **line** not the **field**. You can try:

awk 'BEGIN{OFS="\t";FS="\t"} {split($1,ARR,":"); $1 = ARR[1]; print $0}' inFile.txt > outFile.txt

This will yield:

SNP	Freq.A1	CHR	BP	A1	A2	OR	SE	P
10	0.3519	10	100968448	t	aa	1.0024	0.01	0.812
10	0.4493	10	101574552	a	atg	0.98906	0.0097	0.2585
10	0	10	10222597	a	at	0.9997	0.01	0.9777

The `awk` script is doing a few things:

`BEGIN{OFS="\t";FS="\t"}` is defining the Output Field Separator and Field Separator as tabs.

`split($1,ARR,":")` splits **just the first field** on the colon, and stores each element in the array `ARR`

`$1 = ARR[1]` redefines the first field as the first array element (just the chromosome in this case).

`print $0` prints the entire line (with the newly defined first field).
Question: Hi,

I am currently curating a list of prognostic genes and gene variants responsible for prostate and/or breast cancer. I am currently concentrating only upon germline variants, but am also collecting the somatic variants as when I run into them. Is there a database which already has such a list of genes and also the papers relating to them? The information which I am looking for and curating are the genes, gene locations, mutation type and the link to the paper. Is there any publicly available database which lists all this information.

Thanks,
Venkat

Answer: Check the BROCA panel for breast cancer:

- https://testguide.labmed.uw.edu/public/view/BROCA

For prostate see this paper:

- https://www.ncbi.nlm.nih.gov/pubmed/30777372
Question: Hello!

I'm a newbie in bioinformatics and with QIIME pipeline.

I received 16S rRNA amplicon sequencing data from Illumina MiSeq.
16S rRNA amplicons were sequenced from both ends - paired end sequencing.
Among the different folders that i received, from the core sequencing facility, it was a "raw" folder. This contain individual fastq files for each end (forward and reverse) and for each sample (9 in total). So, the fastq files that i received were already demultiplexed, without barcodes, but with primers.

My questions are:

(1) How do i start? By merging the files and then removing the primers or the opposite?

(2) How do i perform any script that requires the "map.file" if i don't know the barcode sequences?

(3) How do i can perform downstream analysis, such as beta and alpha diversity, without map.file?

I started by running the `multiple_join_paired_ends.py` script. Now i'll try merge all the samples (9) in just one fastq file. Then i was thinking to remove both primer sequences, forward and reverse, through `extract_barcode.py` using the following argument `--input_type -barcode_paired_end`.

Please, help me.
This is really frustrating for someone that just started to learn about bioinformatics and pipelines like QIIME without previous experience on that.
Regards,
@renh@

Answer: (1) How do i start? By merging the files and then removing the primers or the opposite?

Ans:
The first few steps would be:

A. Stitching HQ reads using FLASH or BBMerge

FLASH - https://ccb.jhu.edu/software/FLASH/

BBMERGE- http://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbmerge-guide/

B. Removing chimeric DNA

C. Creating a mapping file.

Start here:

a. http://www.wernerlab.org/teaching/qiime [Tutorial with explanation]

b. [Youtube videos - extensive explanation]

https://www.youtube.com/watch?v=6564K4-_DBI

c. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3249058/ [Research paper]

(2) How do i perform any script that requires the "map.file" if i don't know the barcode sequences?

Ans:
See this link: http://qiime.org/documentation/file_formats.html

See similar links on QIIME1 google group forum:

https://groups.google.com/forum/#!topic/qiime-forum/Lo5DqFdfZh4

https://groups.google.com/forum/#!msg/qiime-forum/XDQnB_QPfHI/xOejWhxwiJwJ

Still, confused, post question on QIIME1 google group

Link: https://groups.google.com/forum/#!forum/qiime-forum

(3) How do i can perform downstream analysis, such as beta and alpha diversity, without map.file?

Mapping file will be required for few steps. Hence, learn how to create mapping file for samples without barcodes. See point (2) as mentioned above.

Try to create a file and run the validate_mapping_file.py script. It generates logs and error that are intuitive. Even, sometimes it generated "corrected" file for you.

Go step by step, it's very very easy.
Question: Hi all,

I am looking for a programmatic way to download complete genomes from RefSeq for a specific taxonomy identifier, in my case I am interested in the Lactobacillales order.

* From the [Genome Download FTP](https://www.ncbi.nlm.nih.gov/genome/doc/ftpfaq/#allcomplete) it is not possible to filter by the order, only genus and species.
* I tried to follow the steps described in this [post](https://ncbiinsights.ncbi.nlm.nih.gov/2013/02/19/how-to-download-bacterial-genomes-using-the-entrez-api/). First two steps (esearch and elink) are okay, but then I do not know how to select my genomes of interest according to my criterion (complete and from RefSeq).
* I also tried the [Ebot pipeline genertaor](https://www.ncbi.nlm.nih.gov/Class/PowerTools/eutils/ebot/ebot.cgi), but then again I am not sure about the query qualifiers to apply.
* My best approach so far is to use the [Genome browser](https://www.ncbi.nlm.nih.gov/genome/browse/#). From there, I can apply my filters and then download the selected records as a .csv or .txt. How can I get the same output file without using the web interface?

Thanks in advance for any hints and help!

Answer: Try this, it should work.

1. get taxids of species belonging to an order using [taxonkit list](http://bioinf.shenwei.me/taxonkit/usage/#list) ,
2. and download all refseq complete geomes according to the `species_taxid` refering to this [post](http://blog.shenwei.me/manipulation-on-ncbi-refseq-bacterial-assembly-summary/) (section "Filter by species_taxid") .
Question: As simply as possible please, and if you don't mind giving a quick definition of any specialist terms you use.

Answer: We really tried to write the main text of [the paper][1] such that it would be understood by non-statisticians. that said, I'll try to do it in a few sentences. before further questions on details please at least *try* to read the paper :)

Let's say we want to compare counts between two groups. We build a model for the observed counts. This model has some parameters: (1) a normalization parameter, for differences in library size at least, or it can be extended by other software; (2) a variance parameter, called dispersion; (3) parameters representing the group differences. Fit (1) using the same method from the original DESeq. Fit (2) in two steps: first find the value of the parameter that makes the [likelihood][2] largest, which is called [maximum likelihood][3] estimation. Look at all the values from all of the genes and move these values towards a middle value. [Bayes theorem][4] guides the amount of movement for each gene: if the [information][5] for the gene is low, the value is moved more to the middle, if the information for the gene is high, the value is moved very little. Fit (3) using the same technique as used for (2). The values for (3) are a useful final product, as are sets of genes where the group differences are likely to be above a threshold (zero or otherwise). These sets are defined by their [false discovery rate][6].

[1]: http://genomebiology.com/2014/15/12/550
[2]: http://en.wikipedia.org/wiki/Likelihood_function
[3]: http://en.wikipedia.org/wiki/Maximum_likelihood
[4]: http://en.wikipedia.org/wiki/Bayes%27_theorem
[5]: http://en.wikipedia.org/wiki/Fisher_information
[6]: http://en.wikipedia.org/wiki/False_discovery_rate
Question: Hi all,

I have downloaded 56G data from ICGC and I am having trouble analyzing it on R (Although I have access to a server). I know that ICGC data is stored in GDC. However, to analyze big data on R, I see suggestions like using SQL through dplyr/R.

Do you know any way to do that? And I would like to know how do you work with big data from databases such as ICGC and TCGA when you need to use R.

I downloaded methylation data for 269 donors `meth_array.PACA-AU.tsv.gz` from here using *wget*:

https://dcc.icgc.org/releases/current/Projects/PACA-AU

I want to compare mean methylation level of selected cg probes across patients.

I would appreciate any help!

Answer: We don't need full dataset in R, subset before importing, for example using *data.table::fread* + *bash* (not tested):

# use head meth_array.PACA-AU.tsv to find out which columns we need
# here we need columns 4 and 9
# 4=icgc_sample_id
# 6=probe_id
# 9=methylation_value

# read only columns 4 and 9 for probe cg00000029
x <- fread("grep -E '^cg00000029$' meth_array.PACA-AU.tsv | cut -f4,9")

# then aggregate as usual
x[ , .(myMean = mean(methylation_value)), by = icgc_sample_id]
Question: Dear All,

I am searching for human 28S, 18S, 5.8S and 5S rRNA (RNA28S, RNA18S, RNA5-8S and RNA5S genes) and also 12S and 16S mitochondrial rRNA (MT-RNR1 and MT-RNR2) in human_genes.gtf file from UCSC.

I am unable to find these genes(RNA28S, RNA18S, RNA5-8S, RNA5S, MT-RNR1 and MT-RNR2 genes) in `human_genes.gtf` file. I am looking at `gene_id` or `gene_name` in gtf file.

How do I check these genes in gtf file or in any other resources.

Is it saved in some other name?

Answer: Hi, guys!

I believe that some 45S genes are now present in the hg38 human reference genome, for example:
https://www.ncbi.nlm.nih.gov/gene/106631777

Here is the location in UCSC GenomeBrowser -- chr21:8,433,222-8,446,572

http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr21%3A8433222-8446572

It is annotated in NCBI, but not in GENCODE (as of version 31).
Question: Hi, everyone,

I am using BLAST+ to identify homologous genes. I would like to know the default parameters of BLASTN. When typing `blastn -h` or` blastn -help`, I could not find detailed default values. On blast+_user_manual.pdf downloaded from NCBI, I could not find these default settings. Could you please do me a favor? I am particularly concern about word_size, penalty score settings, repeat filter, etc. THANK YOU very much!

Answer: 1. **`blastn` (megablast) [for highly similar sequences]**

- word_size = 10
- match/mismatch scores = 1, -2
- gap cost = linear
- filter = low complexity regions

2. **`blastn -task blastn` [quite similar sequences]**

- word_size = 11
- match/mismatch scores = 2, -3
- gap costs = existence: 5, extension: 2
- filter = low complexity regions

You can find this information at NCBI BLAST. If you are running local BLAST with `outfmt` 0 or 1, you will also see the parameters at the end of your output file.
Question: I had thought I understood the meaning of the positive and negative normalized enrichment scores (NES) produced by GSEA, but it just occurred to me I don't know what to expect from a gene set that is enriched both at high and low ranks. For example, some gene sets (e.g., GO biological processes) contain genes that are co-regulated under different conditions, but actually move in opposite directions (producing a mix of positive and negative fold-changes), such that these genes fall at the very top and bottom of the rank-ordered list that GSEA is using to compute NES. So, do these high- and low-ranking genes just cancel each other out leading to an NES close to zero or does GSEA just report the larger absolute enrichment between the positive and the negative NES?

Answer: The short answer to your question is that if you have an equal distribution of up- and down-regulated genes within a gene set your NES (Normalized Enrichment Score) should be close to zero.

The medium answer is that the NES is calculated from the ES (Enrichment Score), which is a rolling sum of the score for each gene. Effectively, you have a ranked list of all genes, and as you move from the top to the bottom of all genes you're changing your ES: for each gene **in** your list you're **adding** to the ES and for each gene **not in** your list you're **subtracting** from the ES. The ES is also weighted based on how close the genes are to one another in the list. So if you have a list of 100 genes, and 50 are right at the top of your list (the 50 most up-regulated genes) and 50 are spread throughout the bottom of your list (all downregulated but with a wide range of log2FC values) then you may still end up with a strongly positive ES.

The long answer would be to dig through the GSEA publications and the (huge amount) of math they use to calculate their scores [here][1] and [here][2].

[1]: https://www.pnas.org/content/102/43/15545.abstract
[2]: https://www.nature.com/articles/ng1180
Question: Hi All,

I would like to ask that is there an accurate tool/script for deinterleaving fastq files? I managed a 'bam to fastq' conversion according to http://gatkforums.broadinstitute.org/discussion/2908/howto-revert-a-bam-file-to-fastq-format and I would like to deinterleave the output fastq.

Many thanks!

Answer: Here is lightweight solution with [Pairfq][1]:

curl -sL git.io/pairfq_lite | perl - splitpairs -i interl.fq -f 1.fq -r 2.fq

Where the input is the interleaved fastq and 1.fq and 2.fq are the forward and reverse reads, respectively (you can choose better names though!). This assumes the pairs are labeled properly (e.g., "/1" or " 1").

[1]: https://github.com/sestaton/Pairfq
Question: Summarized version:

I have an old data set from 2008 from a set of HumanCNV370-Quads, and I have downloaded relatively recent a set of extended-VCFs from the Altai Neanderthal and Denisovan genomes. I want to compare the data between the two. I know the genome coordinates for any given base can shift from assembly to assembly, but will the rsID for a given SNP change if and when the coordinate changes?

In-depth version:

I have SNP data from 64 samples at ~330,000 rs ids (I know there is no mt/Y data, I am pretty sure this is all autosomal). The data is from an old set of HumanCNV370-Quads from 2008. I don't have the genomic coordinates.

I have download two sets of VCF files from [the Denisovan 30x][1] and [Altai Neanderthal 50x][2] coverage genomes. These files are in a cumbersome extended VCF format described [here][3] (page 16) and [here][4] (page 14). These files have rsIDs labeled for most sites (of course though not all sites in these genomes have been assigned rsIDs).

I also have Illumina data from 171 samples (the libraries enriched for NRY- and mtDNA), which I am now have in raw, un-rsID-ed, unfiltered VCFs, which I am trying to bring into mix, but I am going to ignore them for now (I have a thread on them [here][5])

For the 330k sites, I have the alleles for the common chimpanzee from the 1000G (phase1_release_v3/20101123) and the dbSNP build 141 for most of the sites from a friend of a friend. The goal is to use f4 statistics to calculate Neanderthal ancestry estimates.

Anyways, to my question, would the rsIDs from the SNP chip still correspond to the rsIDs that I find in extended VCF files? If not, what would I have to do to make them match up?

[1]: http://cdna.eva.mpg.de/denisova/VCF/hg19_1000g/
[2]: http://cdna.eva.mpg.de/neandertal/altai/AltaiNeandertal/VCF/%29
[3]: http://www.sciencemag.org/content/suppl/2012/08/29/science.1224344.DC1/Meyer.SM.pdf
[4]: http://www.nature.com/nature/journal/v505/n7481/extref/nature12886-s1.pdf
[5]: https://www.biostars.org/p/110272/

Answer: rsID's will stay the same from genome build to genome build, even if the genomic coordinates change. The one gotcha to your plan I can think of is if the rsID was removed from dbSNP or merged with another rsID.
Question: [`vcf2maf`](https://github.com/mskcc/vcf2maf) is the standard tool for converting maf and vcf files back and forth. However, I am encountering situations where the mutation notation between maf and vcf formats are not the same.

For example, I have a mutation in maf format that looks like this:

15 66729115 66729129 + In_Frame_Del DEL GGAACCAGATCATAA -

After converting to vcf format, the mutation is now written like this:

15 66729114 . CGGAACCAGATCATAA C

The `POS` chromosome position value has been changed, along with the ref and alt alleles, to match the requirements of the vcf spec.

What I need is a method for converting maf <-> vcf in such a way that I can retain the original coordinates in the output file. Thus, I want to have an output file that lists something like this, where both the old and new mutation notations are recorded;

15 66729114 CGGAACCAGATCATAA C 15 66729115 66729129 GGAACCAGATCATAA -

-----

For context, I have maf files with extra metadata associated for each variant (which isnt preserved when converting to vcf), and I have converted them to vcf format in order to get more required metadata (in this case, output from `bcftools isec` which only takes vcf input and outputs a `sites.txt` file with presence/absence of each mutation in each sample). Upon trying a simple merge between the two datasets based on the Chrom, Pos, Ref, Alt values in each, I find that I am not able to correctly merge these entries where the mutation notation is change during the conversion. So I need an output file that has both notations recorded to assist in backfilling the new metadata from the vcf analysis into the correct maf variant entries.

Answer: vcf2maf adds a MAF column named `vcf_pos` to preserve the original value in the `POS` field of a VCF. Similarly preserved are `ID`, `QUAL`, and `FILTER` as [documented here][1]. maf2vcf similarly restores these fields into a VCF. If this is the extent of metadata you want preserved between VCFs and MAFs, then you're good to go.

vcf2maf has options `--retain-info`, `--retain-fmt`, and `--retain-ann` to create MAF columns for values in the VCF `INFO`, `FORMAT`, and `ANN` fields respectively. But there is no equivalent `--retain-cols` feature in maf2vcf, that would create `INFO` fields containing MAF data. I have not seen situations where this is needed, that couldn't be solved in other ways.

In all honesty, the MAF format should be obsoleted - and secondary analysis pipelines should only use VCFs in all steps until it is time to report events in cBioPortal or a tab-delimited format for R-users. maf2vcf should not exist. :)

[1]: https://github.com/mskcc/vcf2maf/blob/main/docs/vep_maf_readme.txt
Question: Hi,

Any idea if there's a way to get the Instrument type (i.e. NextSeq, HiSeq, etc), from the instrument name field in the Fastq read header? - I've tried looking for a map of sorts in Illumina's documentation, but to no avail.

Thanks,

Answer: Illumina provided the following information regarding FCID mapping to instruments. The indicated characters are the last four of the nine-character FCID (excepting MiSeq). Please note that this is not comprehensive and there may be new additions in the future.

AAXX = Genome Analyzer
BCXX = HiSeq v1.5
ACXX = HiSeq High-Output v3
ANXX = HiSeq High-Output v4
ADXX = HiSeq RR v1
AMXX, BCXX =HiSeq RR v2
ALXX = HiSeqX
BGXX, AGXX = High-Output NextSeq
AFXX = Mid-Output NextSeq
5 letter/number = MiSeq