cannin/biostars_qa huggingface dataset:
Question: This may not be technically related to bioinformatics.

Using bioinformatics method, lots of biomarker (methylation; microRNA; mRNA; etc...) based classifiers (for benign / malignant tumor; tumor early detection; tumor subtype; etc) have been published and claim to have good performance. However, it seems that few of them have been used in clinical setting. I can only guess that all these classifiers have some kind of shared issue that forbid them being used in clinical setting. Can anyone please share some opinions on this matter?

Answer: Hey, thanks for posting. Other opinions than mine are very welcome.

Some issues that I have identified over the years:

1, lack of reproducibility of results
-------------
This is compounded by poor study design, which itself relates to any number of things:

- no statistical power
- imbalanced sample groups
- failure to control for sources of bias, including age, gender, sampling time, ethnicity, *et cetera*
- untested software that contains bugs

I could add here the fact that every instrument that we use has windows of specification in which the instrument is meant to be operated, and they also have error rates. No NGS instrument, for example, can faithfully sequence any sample of DNA - error will always exist.

2, lack of appreciation of biological variability and how to best capture this
------------
This is mainly for expression, ChIP, metabolomic, and proteomic studies

3, slow and costly clinical trials
-----------------
Clinical trials are very costly and take many years to conduct. Most go nowhere.

4, No 'translational' mechanisms in place
----------------------
In most cases, there are no official mechanisms / systems in place such that data from a research setting can be readily used in a clinical environment, or, if one exists, validation obviously has to take place, and this may be governed by national and / or international law. If something new is introduced, I notice that it is usually a 'local' change to the health system, i.e., in a hospital unit / department, and not something that is global.

------------------------------
--------------

On a side note, I have noticed research change even in the time during which I have been in it (or maybe it is that I have changed). I notice researchers exhibiting greater signs of stress and with less focus on the end goal of the very research that they are conducting. Most do not ponder on what their results could do in terms of improving a health service. Many are also focused intensely on publications and winning the next grant to simply stay in the job. This is a vicious cycle that is ultimately eroding quality research.

I should finally add that I believe more tests will be produced in the next few years. In fact, in certain countries (e.g. Brazil) where there is less regulation, many tests are enrolled straight into practice from research settings. This may prove dangerous, though.

Kevin
Question: Given a missense variant ID, is there a way (via a database or even flat file) to request its estimate frequency or some other importance metric?

By importance metric I mean: how common is this mutation in the average global population or how many times it happens for each 1000 individuals etc.

Thanks you in advance!

Answer: Of course, check the FAQs from the 1000 Genomes paper: http://www.1000genomes.org/faq/how-can-i-get-allele-frequency-my-variant

You can also use the --indv option in vcftools to select only the individuals belonging to a single population, to obtain the MAF only in that population.
Question: I'm tying to extract a specific region of a bam-file into a fasta-file (ultimately). All of the methods I've tried so far give me all reads that OVERLAP the desired region, I'm trying to find a way to trim those to only the desired region.

I've tried:

samtools view

samtools view compiled.sorted.bam ConB:2185-2195

intersectBed

intersectBed -b test.bed -abam compiled.sorted.bam -ubam > out.bam

but these will give the entire read that overlaps my desired region, I'm trying to get something that will trim everything to sam/bam file where the 'reads' are 10 nucleotides long. Am I just missing a flag somewhere to limit the returned region?

Answer: I wrote [SAM4WebLogo](https://github.com/lindenb/jvarkit/wiki/SAM4WebLogo) for https://www.biostars.org/p/73021/ and I think it could do what you need

$ java -jar dist/sam4weblogo.jar -r seq1:80-110 sorted.bam 2> /dev/null | head -n 50
>B7_593:4:106:316:452/1
TGTTG--------------------------
>B7_593:4:106:316:452a/1
TGTTG--------------------------
>B7_593:4:106:316:452b/1
TGTTG--------------------------
>B7_589:8:113:968:19/2
TGGGG--------------------------
>B7_589:8:113:968:19a/2
TGGGG--------------------------
>B7_589:8:113:968:19b/2
TGGGG--------------------------
>EAS54_65:3:321:311:983/1
TGTGGG-------------------------
>EAS54_65:3:321:311:983a/1
TGTGGG-------------------------
>EAS54_65:3:321:311:983b/1
TGTGGG-------------------------
>B7_591:6:155:12:674/2
TGTGGGGG-----------------------
>B7_591:6:155:12:674a/2
TGTGGGGG-----------------------
>B7_591:6:155:12:674b/2
TGTGGGGG-----------------------
>EAS219_FC30151:7:51:1429:1043/2
TGTGGGGGGCGCCG-----------------
>EAS219_FC30151:7:51:1429:1043a/2
TGTGGGGGGCGCCG-----------------
>EAS219_FC30151:7:51:1429:1043b/2
TGTGGGGGGCGCCG-----------------
>B7_591:5:42:540:501/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410/1
TGGGGGGGGCGCAGT----------------
>B7_591:5:42:540:501a/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410a/1
TGGGGGGGGCGCAGT----------------
>B7_591:5:42:540:501b/1
TGTGGGGGCCGCAGTG---------------
>EAS192_3:5:223:142:410b/1
TGGGGGGGGCGCAGT----------------
Question: Hi pretty basic question, but I have a data.frame that looks as follows:

chrom start stop frag dev span
1 chr1 3466334 3466335 60 0.09543313 0.94
2 chr1 4970790 4970791 51 0.08546289 1.00
3 chr1 5022767 5022768 27 0.45680640 0.20
4 chr1 5022807 5022808 34 0.35155312 0.20
5 chr1 5022904 5022905 45 0.12787640 0.20
6 chr1 5023063 5023064 29 0.26792518 0.20

I am trying to do a scatter plot at specific genomic locations, just not sure how to limit it to the right chromosome, here is my ggplot2 command:

my_data %>% ggplot() + geom_point(aes(x = start, y = frag)) + scale_y_continuous(limits=c(0,80)) + xlim(109005904, 109017426) + xlab("Genomic Position") + geom_smooth(aes(x = start, y = frag), method="loess", span=.66, color = "red", fill = "black") + theme_bw()

Thanks!

Answer: Maybe:

my_data %>% filter(chrom == "chr1") %>% ggplot() + ...
Question: I am trying to run blastx (version 2.2.28) with the following command line:

blastx -db nr -outfmt 5 -evalue 1e-3 -word_size 3 -show_gis -num_alignments 1 \
-max_hsps 20 -num_threads 5 -out transcripts_refonly.xml -query transcripts_refonly.fa

So when blastx is running, I can see the output file that grows in size, however when blastx finishes, the output file then is empty (0kb).
This happens when I run it with big files like transcprit fasta file (from a RNA-seq project) of ~30Mb (21000 sequences), but if I run just a smaller file (20 sequences) it works good.

I tried also with .txt ouput format and it happens the same, and also run it into the cluster and different servers and still not working.

Does someone knows what I am doing wrong?

Thank you

Answer: I guess that it was something about memory.. still don't know.
My solution: just split my file into several small files (~100 sequences per file) and run multiple blastx at the same time in the cluster.
Question: Dear community,

I am working with a gene family across a wide range of chordate species, and one of the steps was to retrieve potential homologs from Ensembl. Interestingly, there were some of these genes that had their ID deprecated in newer versions of Ensembl. However, annotated sequences from these deprecated genes had reasonably good quality, and they matched the phylogenetic position that I'd expect, in phylogenetic analyses (so I could tentatively assume that the annotation was correct).
I am curious of the reasons why some gene annotations are deprecated from version to version, even though they look like the annotation is correct, and whether it is wise to use such sequences nonetheless.

All best,

Answer: All these genes are slightly different.

The Mouse Lemur genes were all changed we moved to a new genome assembly. [ENSMICG00000011090](http://aug2017.archive.ensembl.org/Microcebus_murinus/Gene/Summary?db=core;g=ENSMICG00000011090;r=KQ056493.1:1017881-1023258;t=ENSMICT00000011083) and [ENSMICG00000011081](http://aug2017.archive.ensembl.org/Microcebus_murinus/Gene/Summary?g=ENSMICG00000011081;r=KQ056493.1:1033459-1039527) both disappeared when we moved from v2 to v3 of the assembly, whereas [ENSMICG00000011093](http://jul2016.archive.ensembl.org/Microcebus_murinus/Gene/Summary?db=core;g=ENSMICG00000011093;r=GeneScaffold_456:993434-999177;t=ENSMICT00000011088) was lost when we moved from v1 to v2. The good news is that the genes for two of them, [*HTR3C*](http://www.ensembl.org/Microcebus_murinus/Gene/Summary?db=core;g=ENSMICG00000049301;r=1:76618788-76624165;t=ENSMICT00000073294) and [*HTR3E*](http://www.ensembl.org/Microcebus_murinus/Gene/Summary?db=core;g=ENSMICG00000049068;r=1:76660281-76677598;t=ENSMICT00000072625) are in the latest version of Ensembl. The contigs underlying both of them have changed significantly, and this is the reason why the identifiers were changed and not mapped between the releases. I cannot find Mouse Lemur *HTR3D* in the current Ensembl database, although I suspect that [ENSMICG00000049384](http://www.ensembl.org/Microcebus_murinus/Gene/Summary?db=core;g=ENSMICG00000049384;r=1:76604424-76611926;t=ENSMICT00000072460) could be the correct gene based on sequence and genomic position – I will see if this ought to be annotated as this gene.

Kangaroo rat [ENSDORG00000003549](http://dec2016.archive.ensembl.org/Dipodomys_ordii/Gene/Summary?db=core;g=ENSDORG00000003549;r=GeneScaffold_6098:3429-11276;t=ENSDORT00000003549) is different. It was not lost when there was a new genome assembly. Looking that the gene on the assembly at that time, it has two exons that overlaps gaps in the assembly, which is probably why it was deprecated. Since it doesn't have a gene name, I [searched for the sequence using BLAST](http://www.ensembl.org/Dipodomys_ordii/Tools/Blast/Ticket?tl=NdAcgUH9C90qs6W6) in the current database, and there are few solid options for the gene it might be now. It has since had a new genome assembly and the gaps repaired, so could be annotated properly.
Question: I know my samples are from C57BL6.

Should I use the [specific for C57BL6][1] or just [generic mus musculus][2] reference transcriptomes/genomes for alignments?

I am worried that the file size for the reference transcriptome for C57BL6 (`Mus_musculus_c57bl6nj.C57BL_6NJ_v1.cdna.all.fa.gz`, 39.7 MB) is 22% smaller than its generic counterpart (`Mus_musculus.GRCm39.cdna.all.fa.gz`, 51.2 MB). Clearly, biologically the C57BL6 transcriptome cannot be 22% smaller than the transcriptome of some other strain, so there is something going on with the lower detalization for C57BL6??

What do I gain, what do I lose if I opt for C57BL6-specific transcriptome/genome?

[1]: https://ftp.ensembl.org/pub/release-109/fasta/mus_musculus_c57bl6nj/
[2]: https://ftp.ensembl.org/pub/release-109/fasta/mus_musculus/

Answer: There's some basic literature on the topic here: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010552 . Slightly improved uniquely mapped reads; though mapping parameters appear to have a stronger impact than choice of reference.

However, for *differential* expression (BL6/J untreaded vs BL6/J treated) the question isn't so much "how do the quantifications change" but "how do the logFCs change" -- and I don't see published results on this. I should imagine that, by aligning to the strain transcription, some genes might get slightly higher coverage to bump them over the soft-filtering threshold; but few (if any) logFC values should alter based on the reference.

It's probably worth doing twice just to put your mind at ease.
Question: HI!
I would like to get an rs number using (any) of the following information that I have in a data spreadsheet
Gene name (e.g. EXOSC3)/ HGVS cDNA (old nomenclature) (e.g. c.92G>C)/ HGVS protein (p.G31A for protein)/ cDNA_reference (e.g. NM_016042.2)
I was thinking that this would be relatively straightforward - i.e. download similar data with rs number - then merge the files - but can't find this! Would be extremely grateful if anyone can help, Thanks so much, Fleur

Answer: You can use the [UCSC Genome Browser][1] to get this information, but the process depends on the specific HGVS terms you have.

The general idea is to split your list into two lists, one with just your different HGVS terms, and one with just your gene names. Then you can use a combination of the [Table Browser][2] and [Data Integrator][3] to get the rsID's (if any) corresponding to your gene names, and the [Variant Annotation Integrator][4] to get the rsID's (if any) corresponding to your HGVS terms.

We can provide more specific steps if you email this question along with some example data to our mailing list genome@soe.ucsc.edu, that way our entire team will see the question and can help with a solution.

Thanks,

ChrisL from the UCSC Genome Browser

[1]: http://genome.ucsc.edu "UCSC Genome Browser"
[2]: http://genome.ucsc.edu/cgi-bin/hgTables "Table Browser"
[3]: http://genome.ucsc.edu/cgi-bin/hgIntegrator "Data Integrator"
[4]: http://genome.ucsc.edu/cgi-bin/hgVai "Variant Annotation Integrator"
Question: I want to have a bigwig file converted into gff after having the genome coordinates and genome annotation files converted from Hg19 to Hg18. One way, would be to convert bigwig into BED, use the liftover tool and then have the BED file converted into gff. Could someone please suggest a method to convert bigwig to BED, nothing available on the net.

Answer: Instead of bigWigToBedGraph, you might use bigWigToWig, instead, as bedgraph is BED-like but isn't exactly BED, which may affect use of UCSC's liftover tool. In any case, both binaries are available from the UCSC executables download page.

Once you have the file in WIG format, you can convert it to a UCSC BED file with the BEDOPS wig2bed conversion utility.

Alternatively, you could take the bedgraph output and convert it to BED with an added GNU awk step:

$ awk '{ \
if ($1 ~ /^chr/) { \
print $1"\t"$2"\t"$3"\tid-"NR"\t"$4; \
} \
}' foo.bedgraph \
> foo.bed

You should be able to safely run it through liftover at that point.
Question: I have a list of genes. and i want to extract the snps and indels from my VCF file (that i generated using GATK pipeline ) from genes coordinates on . The list of genes coordinates:

Gene Name Accession_no. Start_Position End_Position Strand
Rv0194 NC_000962.3 226878	230462 +

I was looking bedtools but it is asking for .bed format of genes nd as well .bed of bam files. how to do it ? or any other options/tools/scripts?

Like i tried tabix:

bgzip ERR038736_UnifiedGenotyper_variants_raw_snp.vcf
tabix ERR038736_UnifiedGenotyper_variants_raw_snp.vcf.gz
tabix ERR038736_UnifiedGenotyper_variants_raw_snp.vcf.gz AL123456.3:226878-230462 > Rv0194

and this gave me the variants like this:

AL123456.3 227098 . T C 6730.77 . AC=2;AF=1.00;AN=2;DP=172;Dels=0.
AL123456.3 228069 . G A 7132.77 . AC=2;AF=1.00;AN=2;BaseQRankSum=-
AL123456.3 228168 . G C 6682.77 . AC=2;AF=1.00;AN=2;DP=171;Dels=0.

But this is not a vcf file and i can only extract it one at a time. I want to extract all variants against a list of coordinates and store it in a vcf output.

Can anyone help me it this?

Answer: Hello,

you can convert this list into a valid `bed`file by:

$ cut -f2-4 genes.txt|tail -n+2 > genes.bed

You can than take this bed file with `tabix`

$ tabix input.vcf.gz -R genes.bed

fin swimmer
Question: What exactly is "unmappable regions"? My understanding from some google searches is that they are some short regions on the gene that are difficult to map. Is this correct? If so, why are there short region and long region, aren't they randomly splitted?

Thank you.

Answer: In the genome, there is a lot of what is called "repetitive DNA", these are sequences that appear many times throughout the genome. For example LINE1 and Alu are two types of repetitive sequences, that make up a large fraction of the human genome. Naturally, repetitive DNA is processed in sequencing assays like WGS and ChIP-seq, but aligners have a hard time figuring out where the read comes from as the sequence could have originated from many different places. The same thing happens when there are paralogous genes with very similar sequences, the aligner can't exactly distinguish where the sequence originated. This is why in short read sequencing, a lot of reads are discarded from the analysis as we don't know the true genomic origin of those reads. Long read sequencing mostly avoids this problem.
Question: Apologies if this has already been answered, but if I want to update samtools 1.3 to the latest release, what's the best way to go about it? Uninstall the older version first and then do a fresh install, or is there a update command that automatically does that from the terminal?

Answer: Remove the old version and install the new one (or just install the new one over top). But really, as Wouter mentioned, conda is your friend: `conda install -c bioconda samtools`. I recommend making a separate conda environment for each software version, since then you don't have to worry something breaking because you've accidentally up/down-graded some random dependency of another tool.
Question: Is it possible to use `console.log` in CWL javascript expressions? Also I was wondering about `debugger` statements or stepping through the code? Or what is the recommended way to go about debugging with CWL javascript?

Answer: Rabix composer tool editor (http://rabix.io/) has the JS editor integrated. You can simulate tool inputs and it shows what is returned with the expression.
Question: Hi Biostars!,

I am looking to calculate Aneuploidy across a number of tumor samples downloaded from TCGA and ICGC and I had a couple questions regarding the provided copy number data, specifically the segment mean column (which I think is the most relevant data for my cause).

From what I understand from the documentation and previous questions on biostars the values are in Log2 - from this can I read the following into the data?

- A value of 0 should represent CN2
- If so CN2 represents diploidy and variations represent aberrations
- These changes in copy number are a result of the gain or loss of segments of chromosome

I am rather new at CNV work so please let me know if I am completely missing something or barking up the wrong data column.

Thanks for your time, have a great weekend!

Graeme

Answer: You have it about right. Variants represent aberrations, yes. A value of 0 is probably close to CN2 state, but without some detail on the methods, it isn't really easy to tell. Some tumors may have total copy number significantly less than or more than 2, in which case the 0 may not represent CN2 state.
Question: Hi everyone, I was wondering if I can get only the longest trancript in the corresponding column of the Annovar annotation? -onetranscript argument chooses randomly but not the longest one.

So could you please help me? Thanks in advance

Answer: The answer from the author of ANNOVAR is this:

> There has never been a consensus in the field which transcript should
> be used to represent a gene when multiple transcripts are available.
> The most popular approach is to use the longest transcript nowadays.
> However, in the medical genetics field, for certain specific diseases
> and specific genes, there are 'canonical' transcripts that everybody
> uses by default for historical reasons, and you will need to manually
> select this canonical transcript from ANNOVAR output file to
> communicate with the rest of the field.

[source: http://annovar.openbioinformatics.org/en/latest/misc/faq/]

In a way, he is correct, and I feel that the field should start to embrace (and report) multiple transcript isoforms more and more, even with the increased data load. There is too much reporting of variants on isoforms that may have minimal relevance in the tissue of study. Also, for many well-studied genes, like *BRCA1*, we have identified >10 isoforms; whilst, for other less-studies genes, we don't yet understand the alternate splicing patterns of the gene.

Note that VEP does allow you to output the canonical isoform, but to Ensembl the canonical is always the isoform with the longest CCDS: https://www.ensembl.org/Help/Glossary?id=346

On the last point, researchers even disagree about what canonical means. For some it is the highest expressed isoform in the tissue being studied, which may not necessarily be the longest. At least Ensembl's definition is broad-sweeping and covers all tissues.

Kevin
Question: Hi,

I am currently trying to filter a multi-sample vcf file obtained with samtools mpileup. I would like to keep the positions with all samples depths > 4.

I found a way with SnpSift filter, but I have to specify each sample one by one:

cat test.vcf | java -jar SnpSift.jar filter "(GEN[0].DP>4) & (GEN[1].DP>4) & ... &(GEN[66].DP>4)" > filtered.vcf

(I've got 66 samples ... )

Does anyone know how to apply this filtering option to all samples without specifying each samples one by one? (GEN[*] does not seem to work)

Moreover I would like to apply others genotype calling filters on subsets of my samples, is there a way to apply a filter to a specific subset?

Thanks

Answer: How about generating the script?

cat test.vcf | java -jar SnpSift.jar filter `seq 1 66 | awk '{ printf("%s (GEN[%d].DP>4) ",(NR==1?"":" & "), $1);}'`

See also my tool : https://github.com/lindenb/jvarkit/wiki/VCFFilterJS to filter with javascript.
Question: Could you explain the difference between `STAR`, `KALLISTO`, `SALMON` etc. to experimental Biologist/non-bioinformatician.

If possible, the pros and cons of each pipeline.

**Edit below**

***I ask this because three of my colleagues use this 3 difference tools for RNASeq. Basically to answer the same type of biological questions.***

Answer: Regarding the pros and cons of the programs, you might find this article useful:
[https://www.nature.com/articles/s41598-017-01617-3][1]
It finds that kallisto and Salmon produce near identical results, and that STAR (with HTseq for producing gene counts) is less accurate (due to some of the reasons explained above). In terms of speed/memory requirements, the difference between programs is substantial. In recent benchmarking of kallisto vs. STAR on workflows for single-cell RNA-seq
[https://www.biorxiv.org/content/10.1101/673285v2][2]
we found that kallisto was 2.6 times faster than STAR. More importantly, kallisto used much less memory, in some cases 15x less RAM than STAR. This makes it possible to run kallisto on a laptop rather than a server, and facilitates reproducible workflows.

[1]: https://www.nature.com/articles/s41598-017-01617-3
[2]: https://www.biorxiv.org/content/10.1101/673285v2
Question: I have pair end Fastq files from illumina (**Fast1.fastq** and **Fast2.fastq**). Now, how should I proceed further with de-novo genome assembly? I'am completely new to this field of assembling genome sequences? Which open-end tools can be used?

Please help me with the tools and tutorials which can be used for de-novo assembly of sequences.

Answer: If you download the BBMap package, there is a (hopefully) helpful guide in bbmap/docs/guides/PreprocessingGuide.txt

For isolate bacterial assembly with SPAdes, I recommend:

1) Adapter-trimming (you can do quality-trimming at the same time; I suggest a low cutoff, such as Q10)

2) Artificial contaminant filtering

3) Human contaminant removal

4) Error correction

5) Paired-read merging

...then assemble with SPAdes, using both the merged and unmerged reads.