cannin/biostars_qa huggingface dataset:
Question: Hi,

First, thank you for developing the excellent looking GATB library. I'm excited to use it to test out some new ideas, but I had a few basic questions about using the de Bruijn graph. First, I was curious how I retrieve a de Bruijn graph node given the k-mer to which it corresponds. I looked through the documentation, but I couldn't really find any such ability. Does the de Brujin graph in GATB support this (looking up a node by its k-mer?). Second, the release notes for the new version mention that it is now possible to store arbitrary information at node vertices, but I couldn't seem to dig up how to do this in the documentation. Are there any pointers to this?

Answer: Hi Rob,

Thank you for your interest in GATB. Concerning your questions:

1. Our implementation of the de Bruijn graph implies that you normally can't get a Node instance of the graph from an arbitrary kmer value; in fact, you can (a) iterate all the nodes of a given graph (from the information stored in a HDF5 file), and (b) get the neighbors of a given node N from a Bloom filter based structure in memory (but you have to be sure that node N is in the graph, so it should have been got from the all nodes iteration).

The only possibility is when ***you are sure*** that some kmer value corresponds to a node of your graph; then you can use the `Graph::buildNode` method that will return a Node instance from a kmer as an ASCII string (see [here][1]). It is also possible to get a Node instance from the integer value of a kmer; there is no snippet showing this but I can provide some example if you are interested.

If you **really** want to know if a kmer belongs to a graph, you must make yourself a lookup in the kmers stored in the HDF5 file. Although it looks tedious and time consuming, it must be noticed that kmers in the HDF5 files are stored in partitions (so one can look directly in the good partition from the query kmer) and that kmers of a partition are sorted which makes the lookup easy. I think I should add some snippet showing how to do this because some people could really need it.

2. You can indeed tag information for any node now; we use the very nice [EMPHF][2] library for doing it. The idea is to get a unique integer in interval [0,|G|-1] (where |G| is the number of nodes of the graph) from a kmer value of a node; then you can allocate an array of size |G| of any type of your choice and use the EMPH hash function to get the index of the array to be used. Right now, the documentation may lack some examples; you can have a look at the `Graph::queryAbundance` method (file `Graph.cpp`) that shows how to get the coverage of an arbitrary node. Again, I think we should add some snippet about this topic. Note: if you want get the EMPH hash function during the graph building, you have to use `-mphf emphf` (see [here][1])

[1]: http://gatb-core.gforge.inria.fr/snippets_graph.html
[2]: https://github.com/ot/emphf
Question: Hi,

When executing DiscSNP, i got an exeption at the KISSREADS MODULE phase:

#############################################################
#################### KISSREADS MODULE #######################
#############################################################
/home/amelmek/DiscoSnp/build/bin/kissreads2 -predictions discoRes_k_31_c_auto_D_100_P_1_b_0.fa -reads SRR1576548_1.fastq.gz -co discoRes_k_31_c_auto_D_100_P_1_b_0_coherent -unco discoRes_k_31_c_auto_D_100_P_1_b_0_uncoherent -k 31 -size_seeds 26 -index_stride 6 -hamming 1 -genotype -coverage_file discoRes_k_31_c_auto_cov.h5

EXCEPTION: Unable to open bank './ yuw'

there was a problem with kissreads2

**HELP!!!**

Thanks

Amel

Answer: Hi Amel, hi all,

Following an email discussion we found out that discoSnp was called using the following command:

`~/DiscoSnp/run_discoSnp++.sh -r filename1.fastq.gz filename2.fastq.gz -out output_SNP`

As expressed in the doc (*doc* directory) the read set(s) to use should be defined using a *file of file* descriptor, *e.g.*
`~/DiscoSnp/run_discoSnp++.sh -r fof.txt -out output_SNP`

Here is as short copy of the documentation:

- **Case1**: I've a unique read set composed of a unique read file (reads.fq.gz).

*fof.txt* :

> reads.fq.gz

- **Case2**: I've a unique read set composed of a couple of read files (reads_R1.fq.gz and reads_R2.fq.gz). This may be the case in case of pair end sequencing.
*fof.txt* :

> fof_reads.txt

and *fof_reads.txt*:

> reads_R1.fq.gz

> reads_R2.fq.gz

- **Case3** : I've two read sets each composed of a unique read file: reads1.fq.gz and reads2.fq.gz:

*fof.txt*:
> reads_R1.fq.gz

> reads_R2.fq.gz

- **Case4**: I've two read sets each composed two read files: reads1_R1.fq.gz and reads1_R2.fq.gz and reads2_R1.fq.gz and reads2_R2.fq.gz:

*fof.txt*:
>fof_reads1.txt

>fof_reads2.txt

and *fof_reads1.txt*:
> reads1_R1.fq.gz

> reads1_R2.fq.gz

and *fof_reads2.txt*:
> reads2_R1.fq.gz

> reads2_R2.fq.gz

- And so on...
Question: I have genomic positions in a tab delimited file as follows:

Chr Start End
2 10262865 11801950
2 9637403 10927601
2 25141434 27157396
2 31181368 38044662
2 127808684 129276101
2 236957807 238896574
2 237172736 238896574

**First I want to find overlapping segments within these genomic segments in the same file (edit)**. Then I wish to plot these such that the segments do not overlap on the plot. Is there a tool that can aid in::
1 Finding overlapping regions within the same file
2. Visualising the segments in a non-overlapping manner?

My attempts at plotting this in R results in the segments occupying the same position

Answer: If you want to plot them, you can use [karyoploteR](http://bioconductor.org/packages/karyoploteR/). The function [kpPlotRegions](https://bernatgel.github.io/karyoploter_tutorial//Tutorial/PlotRegions/PlotRegions.html) plot the regions non-overlapping by default and will create the plot I think you asking for.

The first two lines read the data file into a GenomicRanges object and change the chromosome names to the ones used by the default genome (UCSC hg19). [plotKaryotype](https://bernatgel.github.io/karyoploter_tutorial//Tutorial/CreateIdeogram/CreateIdeogram.html) creates the plot and [kpPlotRegions](https://bernatgel.github.io/karyoploter_tutorial//Tutorial/PlotRegions/PlotRegions.html) adds the region to the plot. At the [karyoploteR tutorial](https://bernatgel.github.io/karyoploter_tutorial/) you can find more information on how to add more data and customize your plots.

library(karyoploteR)

regions <- toGRanges("regions.txt")
seqlevelsStyle(regions) <- "UCSC"

kp <- plotKaryotype(chromosomes="chr2")
kpPlotRegions(kp, data=regions)

![enter image description here][1]

[1]: https://raw.githubusercontent.com/bernatgel/karyoploteR_questions/master/biostars_367522/plot.png
Question: I'm analyzing structural variants in the 1000 Genomes VCF files and I have a question about the start position. Let's say I have the following CNV:

Pos: 25688926
Chromosome: 1
ID: esv3585526
Ref: G
Alt:
Info: ...END=25700415;SVTYPE=CNV;VT=SV..

What is the start point of this duplication? Does it start at 25688926 and include G (Ref), or does is start after G at 25688927?

Thanks.

Answer: Structural variation annotation is different in VCF. Please refer to the guidelines [VCF format v4.2][1].

For your example:

1	25688926	DUP_gs_CNV_1_25688926_25700415	G	,	.	PASS	SVTYPE=CNV;END=25700415;CS=DUP_gs;AC=101,12;AF=0.02016773,0.00239617;NS=2504;AN=5008;EAS_AF=0.0169,0.004;EUR_AF=0.007,0.001;AFR_AF=0.0492,0.0015;AMR_AF=0.013,0.0014;SAS_AF=0.0031,0.0041	GT

- Chromosome: 1
- Start: 25688926
- End: 25700415
- Type: multiallelic CNV
- Reference Allele: Copy Number 1
- Alternate Alleles: Copy Number 0, Copy Number 2

On genotypes:

0|0: copy number 2
0|1: copy number 1 (1 + 0)
0|2: copy number 3 (1+ 2)
1|2: copy number 2 (0 + 2)
2|2: copy number 4 (2 + 2)

If you want to parse a vcf in perl

#!/usr/bin/perl
open IN, $ARGV[0];
while(){
next if ($_ =~ /^\#/);
chomp $_;
my @r = split /\t/, $_;
my @info = split /\;/, $r[7];
undef my $e;
foreach(@info) { if ($_ =~ /END=/) { $e=$_; $e =~ s/END=//; }}
print $r[0],"\t",$r[1],"\t",$e,"\n";
}
close IN;

In python

#!/usr/env python
import sys
with open(sys.argv[1],'r') as f:
for l in f:
if l.startswith('#'): continue
r = l.rstrip('\n').split('\t')
e=''
for x in r[7].split(';'):
if 'END=' in x: e=x.replace('END=','')
print r[0],r[1],e

[1]: https://samtools.github.io/hts-specs/VCFv4.2.pdf
Question: I downloaded clinical data of BLCA and BRCA cancers from [firebrowse][1] which provides an easy way to access TCGA project's data. Using sample barcodes with sample type code "06"(according to [this link][2]), I could not find any metastatic samples. I also tried to use [GDC Data Portal][3] to find metastatic samples but I failed there too!
Is there any metastatic sample in TCGA project? If yes which cancer types contain more metastatic samples?
Should I use the sample barcode to find metastatic samples?

**Note**: I think this links states the opposite of what I observed!
[TCGA Tissue Sample Requirements: High Quality Requirements Yield High Quality Data][4]
It states that:

> **Sample from primary tumor was necessary**:
A primary tumor is the tumor at the initial site of cancer, not where the cancer may have spread or metastasized, called the secondary tumor. TCGA only studied secondary tumors if they matched to a primary as part of a series (germline, primary and metastasis), or in tumor types for which the primary tumor is rarely diagnosed (such as in melanoma).

[1]: http://firebrowse.org/
[2]: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
[3]: https://portal.gdc.cancer.gov/
[4]: https://cancergenome.nih.gov/cancersselected/biospeccriteria

Answer: I am fairly sure that I have seen samples with code 06 - in fact, almost certain. The number of these samples will differ across cancer and data-type, though. I don't know how samples were selected for profiling - there does not appear to have been any formal agreement among the various centers partaking in the project in this regard.

Note that patients who had a metastatic site are definitively part of TCGA. One can infer this from the clinical data that is held on the GDC Legacy: BRCA

Sample:

barcode	new_tumor_event	new_neoplasm_event_type	new_neoplasm_event_occurrence_anatomic_site
TCGA-A2-A0T7	NO	[Not Available]	[Not Available]
TCGA-A2-A0YC	YES	Locoregional Disease	[Not Available]
TCGA-A2-A0YC	YES	Distant Metastasis	Bone
TCGA-A2-A0YG	NO	[Not Available]	[Not Available]
TCGA-A2-A0YH	NO	[Not Available]	[Not Available]
TCGA-A2-A0YI	NO	[Not Available]	[Not Available]
TCGA-A2-A0YJ	YES	Distant Metastasis	Other, specify

To select out these files, search for your cancer of interest (via the checkboxes at left on the GDC Legacy), and select the *Clinical* data and *Biotab* format.
Question: Hi there,

Crosspost to the google group for ETE3:

I'm trying to install etetoolkit (http://etetoolkit.org/) via conda on MacOS and keep getting the same error. Google has been of very little help unfortunately. I've also removed and reinstalled miniconda.
As suggested on the website, I ran the following:

conda create -n ete3 python=3
conda activate ete3
conda install -c etetoolkit ete3 ete_toolchain

Then I get the following error indicating it cannot find python (weird):

LinkError: pre-link script failed for package etetoolkit::ete3-3.0.0b35-py_0
location of failed script: /PATH/TO/miniconda3/pkgs/ete3-3.0.0b35-py_0/bin/.ete3-pre-link.sh

==> script messages <==

==> script output <==
stdout:
stderr: + /PATH/TO/miniconda3/envs/ete3/bin/python /PATH/TO/miniconda3/pkgs/ete3-3.0.0b35-py_0/link.py

/PATH/TO/miniconda3/pkgs/ete3-3.0.0b35-py_0/bin/.ete3-pre-link.sh: line 2: /PATH/TO/miniconda3/envs/ete3/bin/python: No such file or directory
return code: 127

I've tried putting python in the bin, but that did not solve the problem.

Does anyone have any ideas on how to solve this? I've been using this tool for years with no problem, but it seems like something has gone awry with the latest versions.

Thanks! Please let me know if I can provide more info.
Courtney

EDIT -- SOLVED

Solved!! Try to install the two packages separately, this works!

conda install -c etetoolkit ete3
conda install -c etetoolkit ete_toolchain

Big thanks to the ETE3 team:

https://groups.google.com/g/etetoolkit/c/zH3RLLaSq0Y/m/cLKqtT0DAgAJ

Answer: Solved!! Try to install the two packages separately, this works!

conda install -c etetoolkit ete3
conda install -c etetoolkit ete_toolchain

Big thanks to the ETE3 team:

https://groups.google.com/g/etetoolkit/c/zH3RLLaSq0Y/m/cLKqtT0DAgAJ
Question: I've been calculating differential expression for two separate groups both of which have baseline and treatment expression measurements. I have looked at differential expression between the baseline and treatment and now have 2 lists of differentially expressed genes. Differentially expressed genes were identified with edgeR with appropriate correction for FDR.

It's an interesting question to compare the lists of differentially expressed genes and I can identify genes significantly up- and down-regulated in both groups. I can also show that the overlap in these lists is significant using Fisher's exact test.

However, when I look at genes that are significantly differentially expressed in only one group do I need to do any further analysis or test to ensure the difference is "real"? For example, one gene could be significantly DE in one group close to the cutoff for significant (e.g. P = 0.049) but just misses the cutoff in the second group (e.g. P = 0.051). I've been unable to find anything about this online so far.

Answer: I'd encourage you to use GSEA or a similar rank-based method for comparisons rather than choosing some p-value cut-off and comparing lists.
Question: I am using bowtie2, and due to the size of the SAM output, I would like to entirely avoid the SAM file, and pipe the results into samtools to convert the output to a BAM file. I am getting a little bit hung up on how to accomplish this.

I have tried a few different options:

The approach here: [http://slhogle.github.io/2014/bowtie-and-samtools/][1]

and here: [http://www.metagenomics.wiki/tools/samtools/breadth-of-coverage][2]

are both giving me an empty BAM file.

I have also tried:

bowtie2 -x [index] -U [reads] -S ${TMPDIR}/temp.sam | samtools view -bS > out.bam

To see if I would get anything different using the temp directory to store the SAM temporarily.

Edit in response to a comment to add the commands I have tried:

bowtie2 -x [index] -U [reads] -S - | samtools view -bS > out.bam

And

bowtie2 -x [index] -U [reads] -S - | samtools view -bS - > out.bam

I'm sure this has been done before, is there an obvious way to accomplish this task?

Thanks.

[1]: http://slhogle.github.io/2014/bowtie-and-samtools/
[2]: http://www.metagenomics.wiki/tools/samtools/breadth-of-coverage

Answer: You shouldn't need to use the -S command for the output, just omit this arguement

```
bowtie2 -x genomeindex -U out.fq|samtools view -bS - > out.bam
```

Versions

```
% bowtie2 --version
/usr/bin/bowtie2-align-s version 2.3.4.2
64-bit
% samtools --version
samtools 1.7
Using htslib 1.7-2
Copyright (C) 2018 Genome Research Ltd.
```
Question: I want to receive the protein sequence of the following gene IDs and output a fasta file with the sequences with its identifier.

```
handle = Entrez.esearch(db="gene",
term="primate[Orgn] AND TNF[Gene Name]",
idtype="acc",
retmax='50',
)
record = Entrez.read(handle)
idlist = record['IdList']
print(idlist)
```

But I am not sure where to go from here. Any help would be appreciated.

Answer: Using command line [EntrezDirect][1] (truncated for space) :

$ esearch -db gene -query "primate [orgn] AND TNF [gene]" | elink -target protein | efetch -format fasta > tnf.fa
$ more tnf.fa
>sp|Q19LH4.1|TNFA_CALJA RecName: Full=Tumor necrosis factor; AltName: Full=Cachectin; AltName: Full=TNF-alpha; AltName: Full=Tumor necrosis factor ligand superfamily member 2; Short=TNF-a; Contains: RecName: Full=Tumor necrosis factor, membrane form; AltName: Full=N-terminal fragment; Short=NTF; Contains: RecName: Full=Intracellular domain 1; Short=ICD1; Contains: RecName: Full=Intracellular domain 2; Short=ICD2; Contains: RecName: Full=C-domain 1; Contains: RecName: Full=C-domain 2; Contains: RecName: Full=Tumor necrosis factor, soluble form; Flags: Precursor
MSTETMIQDVELAEEALPKTRGPQGSKRRLFLSLFSFLLVAGATALFCLLHFGVIGPQKDELSKDFSLIS
PLALAVRSSSRIPSDKPVAHVVANPQAEGQLQWLNRRANALLANGVELRDNQLVVPSEGLYLVYSQVLFK
GQGCPSNFMLLTHSISRIAVSYQAKVNLLSAIKSPCQRETPQGAKTNPWYEPIYLGGVFQLEKGDRLSAE
INLPDYLDLAESGQVYFGIIGL
>sp|P48094.1|TNFA_MACMU RecName: Full=Tumor necrosis factor; AltName: Full=Cachectin; AltName: Full=TNF-alpha; AltName: Full=Tumor necrosis factor ligand superfamily member 2; Short=TNF-a; Contains: RecName: Full=Tumor necrosis factor, membrane form; AltName: Full=N-terminal fragment; Short=NTF; Contains: RecName: Full=Intracellular domain 1; Short=ICD1; Contains: RecName: Full=Intracellular domain 2; Short=ICD2; Contains: RecName: Full=C-domain 1; Contains: RecName: Full=C-domain 2; Contains: RecName: Full=Tumor necrosis factor, soluble form; Flags: Precursor
MSTESMIRDVELAEEALPRKTAGPQGSRRCWFLSLFSFLLVAGATTLFCLLHFGVIGPQREEFPKDPSLI
SPLAQAVRSSSRTPSDKPVAHVVANPQAEGQLQWLNRRANALLANGVELTDNQLVVPSEGLYLIYSQVLF
KGQGCPSNHVLLTHTISRIAVSYQTKVNLLSAIKSPCQRETPEGAEAKPWYEPIYLGGVFQLEKGDRLSA

If you want to save individual sequence in a separate file then use:

$ esearch -db gene -query "primate [orgn] AND TNF [gene]" | elink -target protein | efetch -format acc | xargs -n 1 sh -c 'efetch -db protein -id "$0" -format fasta > "$0".fa'

[1]: http://bit.ly/entrez-direct
Question: I'm using IGV and I grouped the alignments by supplementary tags
Colored reads on the top panel are supplementary
![Colored reads on the top panel are supplementary][1]

How do I find the primary alignment of this read?

[1]: http://i.imgur.com/qi6Ba1u.png

Answer: Hi, Zev is correct that IGV wasn't originally designed for SVs, however improvement in this area is a major thrust of development now. Specifically supplementary tags are recognized, and you can link supplementary alignments by right clicking the alignment menu and selecting "Link supplementary reads". If all reads are in view they will be linked similar to paired alignments by a line, or if they overlap you will see this with alpha transparency. Also, linked alignments that are inconsistent with respect to strand will be highlighted with color.

There are also features for 10X data, specifically if 10X barcodes are found in the data you will see 2 menu items, to link or by molecular ID. If phasing information is present they linked alignments will be grouped and colored by phase. This view is somewhat similar, by design, to the Loupe linked reads view.

This effort is still early, and I welcome feedback. The best place to do that is on our git forum, https://github.com/igvteam/igv/issues. There are topics open now specifically for PacBio and 10X, https://github.com/igvteam/igv/issues/277 and https://github.com/igvteam/igv/issues/275. I sometimes don't see postings here unless they are specifically called to my attention.

You can access the new features in the snapshot build, which requires Java 8, at http://software.broadinstitute.org/software/igv/download_snapshot
Question: Hi All,

I want to convert a matrix in numeric format, 0 for homozygous reference allele, 1 for heterozygous and 2 for homozygous alternative to a genotype matrix, if it is possible to a vcf file. Do you know any tool to do it?

Best Wishes,

Eric

Answer: If your matrix is variant-major (i.e. each line has one variant, each column has one sample), it shouldn't take too much work to manipulate it into a form that plink 2.0 --import-dosage (see http://zzz.bwh.harvard.edu/plink/dosage.shtml for a description of the expected format, and http://www.cog-genomics.org/plink/2.0/input#import_dosage for --import-dosage syntax) can read; adding "--export vcf" to the same command line will then give you your VCF.

If it's sample-major, you'll probably need to transpose it first.
Question: I have a vcf file from a single cell experiment. I got the basic statistics using vcf-stat utility in vcf-tools. I am interested in the distribution of indels. How can i get a barplot of genome-wide indel distribution (number of indels per chromosome) as well as a plot of number of insertions/deletions versus size of insertion/deletion? Is there a tool that I can use ?

Answer: The VCF will have all of this information. One way to summarize all of your data would be to use ```grep``` to thin your VCF down to just indels, then read the (tab delimited) result into R and plot whatever you want. Note that, per record, you'll have a chromosome, start position, and stop position, which should give you everything you'll need. Other tools might be able to handle this, but it's just basic text parsing to get what you want.
Question: I want to grab select packages from the Bioconductor repository index. As I understand it, there are one of two ways to do this:

Use some argument with the following repository url:

repos = biocinstallRepos()[1]

Use some extra argument to the makeDepGraph function below:

z = makeDepGraph(repos)

How can I specify - by name - for which packages I want to create a dependency graph?

Also, using Rgraphviz, how can I add the name of each package as a label?

TIA.

Answer: For your first question, take a look at the subGraph function from the graph package. For your second, please take a look at the Rgraphviz vignette.
Question: I am a clinician, and have no particular expertise in genomics. I'm confused by much of what I see on dbSNP. Any help would be greatly appreciated! Thank you in advance.

I need some help understanding HGVS Names!

I see SNP's referred to with conventions like **A66G** for MTRR ([rs1801394][1]) and **A1298C** for MTHFR ([rs1801131][2]).

When I look under **HGVS names** for rs1801394 I see: NM_002454.2:c.66A>G, which l assume equals "A66G", but there are many other "names" there. What are these names? Are they all equivalent?

When I look under **HGVS names** for rs1801131 I do not see anything resembling "A1298C". Does this mean I have the wrong rs# for this SNP?

According to [SNPedia][3], it's the correct rs#: "rs1801131 is a SNP in the MTHFR gene, representing an A>C mutation at mRNA position 1298".

[1]: http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1801394
[2]: http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1801131
[3]: http://www.snpedia.com/index.php/Rs1801131

Answer: the 1298 in A1298C must be relative to the mRNA, which includes the 5' UTR, or maybe just based an earlier prediction as to where the CDS began. I would imagine there is some historical reason it is still known by that name.

Here is one of many papers where they discuss c.665C>T and c.1286A>C: http://www.nature.com/gim/journal/v15/n2/full/gim2012165a.html

HGVS coding nomenclature starts from the CDS start, and is transcript specific, hence NM_005957.4:c.1286A>C is the HGVS identifier for this mutation. As long as you include the c-dot people will know what you mean.
Question: Hi,
I am working on variant calling with a fungal genome. When I am trying to reorder my .bam file (sorted) using picard ReorderSam (picard v2.9.0) to be used further in GATK realignment & BQSR, it throws this error:

Exception in thread "main" picard.PicardException: Discordant contig lengths: read scaffold1 LN=4751343, ref scaffold1 LN=4340956

Command I used,

java -jar /path/to/picard.jar ReorderSam I=myfile.bam O=reordered.bam R=myreference.fa

I checked my reference.fa file and found that the length of scaffold1 mention was incorrect (4751343), so I corrected it (LN=4340956) and created the .dict file again. But, I am still getting the same error.

Answer: Use option ALLOW_CONTIG_LENGTH_DISCORDANCE=true

> If true, then permits mapping from a read contig to a new reference contig with the same name but a different length. Highly dangerous, only use if you know what you are doing