cannin/biostars_qa huggingface dataset:
Question: Hi everyone,

I have a list of 150 locations (eg: "10:112834:113150") and I would like to retrieve the variants found in these locations,

This is the command I use to get the variants (*fantom_locs* contains the list of locations):

snp_db <- useEnsembl(biomart="snp", dataset="hsapiens_snp")

snps_FANTOM <- getBM(attributes = c('refsnp_id'),
filters = c('chromosomal_region'),
values = fantom_locs,
mart = snp_db)

It works for a few locations but then I get a timeout error, I can implement a for loop to fetch the snps on one location at a time, but it is very slow and I just wondered if there was a better way to do this.

For example is it worth it to download a local version of the mart ?

Answer: Hi Corentin,

Since this query is too large for BioMart, I would suggest using the Ensembl REST API Overlap endpoint, restricting the query using the 'feature=variation' parameter:
https://rest.ensembl.org/documentation/info/overlap_region

e.g
https://rest.ensembl.org/overlap/region/human/7:140424943-140624564?content-type=application/json;feature=variation

Best wishes

Ben
Ensembl Helpdesk
Question: Is there any tools available for smoothing the wig file (fixed or variable) with different span (wig file can have with any spanned data) ? Although, [DANPOS][1] has a function for smoothing the resultant wig file, but it does not have any separate function for smoothing the wig file.

Here I am showing an example of before and after smoothing:
![enter image description here][2]

Green: After smoothing
Black: Before smoothing
Note: both files are different

[1]: https://sites.google.com/site/danposdoc/tutorial/dpos
[2]: https://lh3.googleusercontent.com/-XuWEuiub9v4/WbmA_O3MqdI/AAAAAAAAIbE/NxPh2hzlgpIgbt4fcdYlLa0jOzxnCqoEgCL0BGAYYCw/h286/2017-09-13.png

Answer: Here's a way to do things on the command line, which gives quantitative control over the inputs and outputs.

1) Convert your fixed- or variable-step WIG to BED with [BEDOPS `wig2bed`][1]:

$ wig2bed < signal.wig > signal.bed

2) Fetch chromosomes for your genome of interest using Kent Utilities `fetchChromSizes` and create a sorted BED file with BEDOPS `sort-bed`:

$ fetchChromSizes hg38 | awk -v OFS="\t" '($1!~/_/){ print $1, "0", $2 }' | sort-bed - > hg38.bed

3) Use BEDOPS `bedops` with the `--chop` option to generate bins of the desired width.

For example, we can generate disjoint windows that are 100kbases wide:

$ bedops --chop 100000 hg38.bed > hg38.100k.bins.bed

4) Finally, smooth the signal in the BED file with BEDOPS `bedmap`, using the signal BED file from step 1, the bins made from the chromosomes acquired from step 3, and some range that collects signal upstream and downstream of the window.

For example, for each bin, we can calculate the `mean` signal that lies over the window, as well as 50kbases outside the window:

$ bedmap --echo --mean --range 50000 hg38.100k.bins.bed signal.bed > answer.bed

The `answer.bed` file will contain the coordinates of the bin and the average signal over the bin and 50kbases around its sides.

The smaller the window and the wider the range, the smoother the resulting post-processed signal. Like using coarse sandpaper on a fine hardwood, you could sand away too much. You could use a chromosome as sample input to see how tuning the different parameters changes the smoothing.

[1]: http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/wig2bed.html
Question: Hi all,

I received structural variant calling results for the comparison of two bacterial genomes and I have multiple structural variants named "**Breakends**" and characterized as "SYMBOLIC" type. I've tried to search litterature article about it but it seems that there is not a proper definition.

Could I have some help to understand what a breakend is actually?

Thank you in advance,
Audrey

Answer: The best place to help is probably to look at section 5.4 of the VCF specification https://samtools.github.io/hts-specs/VCFv4.3.pdf

This describes it in detail, and it is pretty accessible to read. That said, breakends can be difficult to properly interpret. It can require putting together multiple records (e.g. multiple lines) of the VCF file to understand what is going on.

For example in VCF, you could use the syntax to represent a inversion (this example is from the spec above)

```
#CHROM POS ID REF ALT QUAL FILTER INFO
2 321682 INV0 T 6 PASS SVTYPE=INV;END=421681

```

or that same inversion could be described the breakends:
```
#CHROM POS ID REF ALT QUAL FILTER INFO
2 321681 bnd W G G]2 : 421681] 6 PASS SVTYPE=BND;MATEID=bnd U;EVENT=INV0
2 321682 bnd V T [2 : 421682[T 6 PASS SVTYPE=BND;MATEID=bnd X;EVENT=INV0
2 421681 bnd U A A]2 : 321681] 6 PASS SVTYPE=BND;MATEID=bnd W;EVENT=INV0
2 421682 bnd X C [2 : 321682[C 6 PASS SVTYPE=BND;MATEID=bnd V;EVENT=INV0
```

You might ask, why would you represent it using breakends if you could have the more simple representation? One reason might be that it might not just be a simple INV, where two points in the genome are flipped around...it may be more complex and lead to something like an "INVDUP"

See https://kero.hgc.jp/examples/CLCL/hg38/index.html for a good example and visualization of an INVDUP

![enter image description here][1]

That above thing would likely be represented as a series of breakends. But you could also ask, well, why not just say that INVDUP example like this?
```
2 321682 INV0 T 6 PASS SVTYPE=INVDUP;END=421681
```

The reason might be that breakends can express more precision than that, giving the precise breakpoints, at the cost of being a bit harder to interpret. We have been working on jbrowse 2 to try to make tools for people to interpret breakends, but it is still challenging.

[1]: /media/images/37a716ad-7723-4c0b-8666-f6d80c3d
Question: As referenced in the pvcf documentation… what are overlapping variants?
https://www.biorxiv.org/content/biorxiv/early/2018/06/11/343970.full.pdf

> **Hybrid allelic representation.** To facilitate
downstream summary statistics without doublecounting, ideal unified sites would be completely nonoverlapping, with mutually-exclusive alleles.

I understand snv, insertion, deletion, copy number… but what does an overlap look like? I can’t find any descriptive information aside from mentions that “variants can overlap.”

Is this when a single mutation is represented multiple ways by different ordering of alleles?

Answer: here are two variants from gnomad:

1	13459	rs1315857414	CAGA	C
1	13462	rs1441058751	AAGT	A

as you can see, the first variant CAGA at 13459 overlaps with th position 13462
Question: I have to find out the size of the protein sequence, but even using the codes below, I couldn't. This first code was to find how many proteins there are in total and to find the size of the sequences.
The attached image is just to show what I want the code to search for. I don't know what is missing in the code

arq = open("genoma9.faa")
conteudo = arq.read()
print(conteudo)
fh = open("genoma9.faa")
n= 0
for line in fh:
if line.startswith(">"):
n+= 1
print(line)
proteins = line.count(">")
print("Total of Proteins: " + str(proteins))

![enter image description here][1]

Trying to find this middles characters above the >WP:

Example:

>WP_013277001.1 DNA polymerase III subunit beta [Acetohalobium arabaticum]
MQIKIDRKNFYDGIQTVRKAISSKSTLPILSGILIETQEKKLKLVGTDLELGIECRVDANIIKDGAIVLPANHLANIVRE
LPNKELELELKKDNKIEISCGLSQFKIHGSPADEYPLLPEVGSGIEYTLSQEKFQAMINRIKFATSDDESRPFLTGGLLS

[1]: /media/images/21b38c5b-d255-4d16-8f9c-8bb33815

Answer: I think you are reinventing the wheel. There is no need to write separate code for handling biological sequences when it all exists in BioPython and can be accessed in several lines of code. What I show below could be optimized, so it is only for illustration. I suggest you save it into a file `fasta_len_and_number.py` or something like that.

import sys
from Bio import SeqIO

# open the file specified after script name
FastaFile = open(sys.argv[1], 'r')

counter = 0 # initialize sequence counter
for rec in SeqIO.parse(FastaFile, 'fasta'):
counter = counter + 1 # increase sequence counter
name = rec.id # sequence header
seq = rec.seq # protein/DNA sequence
seqLen = len(rec) # determine sequence length
print(seqLen, name) # print the length + header

print('\n A total of %d sequences' % counter)
FastaFile.close()

Running this line:

python fasta_len_and_number.py genoma9.faa

will make hopefully a desired output. On one of my files when I tested, the last 10 lines look like this:

130 2HY5_A
114 2Q68_A
141 4CDL_A
47 6O3S_A
12 4AKT_C
37 4Z80_C
145 1Y23_A
215 6S5A_L

A total of 112217 sequences
Question: Hi,

I am trying to get the distribution of fragment lengths from a BAM file. Any ideas which programs/tools can do that?

Answer: Ok this is template length according to sam format. Column 9 stores this value.

samtools view Your.bam | awk '{print $9}' | sort | uniq -c
Question: I am planning to use [Bioconductor GEOquery][1] to download a couple of micro-array datasets from NCBI GEO.

[1]: http://www.bioconductor.org/packages/release/bioc/html/GEOquery.html

Then I would like to export a subset of the metadata and the expression data to flat files that I can import elsewhere.

What I have so far is:

```
library(GEOquery)
library("R.utils")

geo_id <- "GSE45016"
gse <- getGEO(geo_id,GSEMatrix=FALSE)

#show metadata
Meta(gse)

#show metadata for first sample
GSMList(gse)[[1]]

#select specific field from metadata of first sample
GSMList(gse)[[1]]@header$characteristics_ch1

# Result for sample 1
[1] "tissue: normal prostate (NP) epithelial cells"

GSMList(gse)[[2]]@header$characteristics_ch1

# Result for sample 2
[1] "tissue: prostate cancer cells" "clinical stage: clinical T4N0M1"
[3] "gleason score: GS 9" "psa level: PSA 5477ng/ml"
```

As you can see the number of key value pairs is different for sample 1 and 2. What is would like to have is an array for every key under

@header$characteristics_ch1

and then the value or null (in case the key is missing) for every sample in the GEO dataset" ;

```
key_tissue: normal prostate (NP) epithelial cells\tprostate cancer cells
key_psa_level: null\tPSA 5477ng/ml
```

Other metadata fields like "title" luckily only have a single value beneath it.

```
GSMList(gse)[[1]]@header$title = "Normal prostate"
GSMList(gse)[[2]]@header$title = "High-grade PC1"
```

Also these I would like to have in an array for the key title.

My second question is how to export the expressions data that is stored under every sample. I would like to stream trough all the probes, get the expression values for that probe for each sample and write it to another csv file.

Answer: I think that the way you have chosen to read the GSE data into R has created some confusion for you.

Try this instead (note: formatting was lost here so posted as a Gist):

https://gist.github.com/neilfws/1dff7fceb641dba374db#file-gistfile1-r

As for exporting the expression data:

exp <- exprs(gse)

returns a matrix where the column names are sample names.
Question: I installed `nextflow` with `conda install -c bioconda nextflow` on HPC, which has Singularity3.7.1 installed.

When I `nextflow -Dcapsule.log=verbose` There is an error:
```
CAPSULE: Jar: /lila/home/lid/miniconda3/envs/nf-core/share/nextflow/dist/20.10.0/nextflow-20.10.0-one.jar
CAPSULE: Platform: linux
CAPSULE: Loading caplet: MavenCapsule
CAPSULE: Applying caplet MavenCapsule
CAPSULE EXCEPTION: Could not read pom.xml while processing system property capsule.log: verbose
java.lang.RuntimeException: Could not read pom.xml
at MavenCapsule.createPomReader(MavenCapsule.java:249)
at MavenCapsule.createPomReader(MavenCapsule.java:242)
at MavenCapsule.finalizeCapsule(MavenCapsule.java:70)
at Capsule.finalizeCapsule(Capsule.java:819)
at Capsule.(Capsule.java:744)
at NextflowLoader.(NextflowLoader.java:37)
at jdk.internal.reflect.NativeConstructorAccessorImpl.newInstance0(Native Method)
at jdk.internal.reflect.NativeConstructorAccessorImpl.newInstance(NativeConstructorAccessorImpl.java:62)
at jdk.internal.reflect.DelegatingConstructorAccessorImpl.newInstance(DelegatingConstructorAccessorImpl.java:45)
at java.lang.reflect.Constructor.newInstance(Constructor.java:490)
at Capsule.newCapsule0(Capsule.java:5603)
at Capsule.newCapsule(Capsule.java:5595)
at Capsule.myCapsule(Capsule.java:361)
at Capsule.main0(Capsule.java:387)
at Capsule.main(Capsule.java:379) Caused by: java.util.zip.ZipException: invalid entry CRC (expected 0xc952da5c but got 0x7155fe57)
at java.util.zip.ZipInputStream.readEnd(ZipInputStream.java:408)
at java.util.zip.ZipInputStream.read(ZipInputStream.java:197)
at java.util.jar.JarInputStream.read(JarInputStream.java:207)
at java.util.zip.ZipInputStream.closeEntry(ZipInputStream.java:141)
at java.util.zip.ZipInputStream.getNextEntry(ZipInputStream.java:119)
at java.util.jar.JarInputStream.getNextEntry(JarInputStream.java:142)
at Capsule.getEntry(Capsule.java:3839)
at Capsule.getEntryInputStream(Capsule.java:3835)
at MavenCapsule.createPomReader(MavenCapsule.java:246)
... 14 more
Unable to initialize nextflow environment ```

What should I do to fix this problem?
Thanks

Answer: > : invalid entry CRC (expected 0xc952da5c but got 0x7155fe57)

looks like there is a problem in a zip/jar file which is badly packaged or installed.

You should try by downloading nextflow itself rather than using your conda installation.
Question: Hi,

I have one .vcf file of whole genome sequencing of tumour Vs normal samples of 21 patients.

I need a data from like this as input for a tool for finding driver genes

> head(mutations)
sampleID chr pos ref mut
1 Sample_1 1 871244 G C
2 Sample_1 1 6648841 C G
3 Sample_1 1 17557072 G A
4 Sample_1 1 22838492 G C
5 Sample_1 1 27097733 G A
6 Sample_1 1 27333206 G A

In separated .vcf files for each patient I have start, end, chromosome, ref, and variant allele. However I am sure how to get such data frame from this big vcf

Any help please?

Thank you

Answer: bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%ID]\n' c.vcf

This solve the error

`bcftools query plugin` and `snpsift plugin` in `galaxy` also do that
Question: Hello,

I will analyze soon RNA-Seq samples build using a TruSeq RNA Exome (illumina) from FFPE. I've already a standard RNA-Seq pipleline that works well on fresh tissue ( RiboZero kits). Pipeline is a classic STAR-DESeq2 + additional post-hoc analysis (splicing, fusion, pathways, etc..).

In this context, should I be carefull on some steps using capture-based RNA-Seq samples from FFPE samples ?

thank you

Answer: FFPE samples often suffer from low quality RNA inputs, i.e. lots of degraded RNA, which often leads to pronounced 3' bias and sometimes dramatically reduced library complexity (i.e., only few genes were captured).

You should make sure to check the gene body coverages, exon coverages and library complexity (e.g. using [QoRTs](https://hartleys.github.io/QoRTs/)) because they may differ quite significantly from the fresh samples and may influence the types of comparison you can do.

Some nice assessments of the impact of degradation were done by [Schuierer et al., 2017](https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-017-3827-y) and [Li et al., 2014](https://www.ncbi.nlm.nih.gov/pubmed/25150837).
Question: I am updating the SICER algorithm ( https://github.com/endrebak/epic ) and have added paired end support. The original SICER does not support paired end reads, but for the single end case it reduces each read to a single coordinate, which is the start of the read plus half the fragment size.

I have no paired end ChIP-Seq data myself to try it on, so I might have implemented a too naive paired end mode. It reduces each paired end read to a point by taking the leftmost and rightmost coordinate of (both the starts and ends of) a pair and finding its midpoint.

So for the (fake) read pair

chr7	20246668	20246669	chr7	20246693	20246694	U0	0	+	+

the coordinate is

20246694 + (20246694-20246668)/2 = 20 246 707

It seems like this might lead to a problem though:

If the two mates are very far apart, the midpoint might be in a bad (ie heterochromatic) or uninteresting (ie blacklisted) region.

What is the best way to solve this?

I can think of two solutions, but do not understand all their up- and downsides:

1) Discard read pairs more than say 100 bp apart
2) Treat each mate in a pair as an individual read

I lean towards the first solution since it seems like the paired end libraries my users use contains much much more data than a typical single end library, and doubling that amount seems like it would be a lot of pain (waiting) for little gain (better results).

Answer: Option #2 looks bad to me, exactly for the reasons you pointed out.

Concerning option #1, the cutoff (100 bp in your example, which is very stringent IMO) could be set independantly for each library or each experiment. Indeed, the expected fragment size can vary depending of experimental procedure (sonication, size selection, ...) For instance :

cutoff = median(fragments_size) + 2*sd(fragments_size)

Finally, even if you don't own PE data, you could download some on SRA and test the robustess of your method on various datasets.
Question: Hi,

HTseq need gff file for reads counts from alingment sam file. I have GTF file which I used for read alignment, How can I convert this into gff file.

Thanks

Answer: You can use GTF for HT-seq, you don´t have to convert nothing.
Question: I was trying to map a list of genes to Oryza sativa (rice) KEGG pathways using [pathview][1] R package. But failed to map genes. I followed [this][2] tutorial.

For example, the code below

library("pathview")
geneList<-read.table("genes.txt",header=T,sep="\t")
pathway <- pathview(gene.data = geneList,
pathway.id = "00195",
gene.idtype = "RAPDB",
species = "Oryza sativa japonica"
)

with the input "genes.txt"

Os03t0659266-00
Os04t0473150-00

Should map input genes to a specific pathway. But showing error

Error in pathview(gene.data = geneList, pathway.id = "00195", gene.idtype = "RAPDB", :
No proper gene annotation package available!

However, the online version shows mapping with the same input. So how can I solve this problem?

[1]: https://bioconductor.org/packages/release/bioc/html/pathview.html
[2]: https://bioconductor.org/packages/devel/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

Answer: Your gene IDs (eg Os03t0659266-00) are not Entrez or KEGG gene ID, hence can’t be mapped to KEGG osa pathways. That’s why you saw the warning and plain graph with no data. you can generate simulated osa data with entrez/kegg ID as below, and it will work. So you need to convert your gene IDs to entrez/kegg gene IDs manually.

>gdata.osa=sim.mol.data(mol.type="gene", species ="osa", id.type="kegg", nmol=10000)
>pv.out=pathview(gene.data = gdata.osa, pathway.id = "00195", gene.idtype = "KEGG", species = "Oryza sativa japonica")

FYI, you can check the ID type supported by pathview for your target species (i.e. osa) using korg data in pathview like:
>data(korg)
> ri=korg[,"kegg.code"]=="osa"
> korg[ri,]

ktax.id tax.id kegg.code
"T01015" "39947" "osa"
scientific.name common.name entrez.gnodes
"Oryza sativa japonica" "Japanese rice" "1"
kegg.geneid ncbi.geneid ncbi.proteinid
"4349630" "4349630" "XP_015615283"
uniprot
"A0A0P0XY71"

So supported gene ID for osa include: kegg.geneid, ncbi.geneid (entrez), ncbi.proteinid, uniprot.
The simulated data with Entrez/KEGG IDs look like:
> head(gdata.osa)

4328398 4326389 4338593 107277186 4334918 4335759
-1.4462875 0.3158558 -0.3427475 -1.9313531 0.2428210 -0.3627679
Question: Hi,

I find this tool [https://github.com/kblin/ncbi-genome-download][1] in github for downloading the bacteria genomes.

ncbi-genome-download bacteria

it downloades ---> refseq/bacteria/GCF_940077525.1

from GCF_ * how do I get fasta files ?

Then I tried this, cmd, it unables to download the fasta files.

ncbi-genome-download --formats fasta bacteria --parallel 16

WARNING: Skipping entry, as it has no ftp directory listed: 'GCF_023646435.1'

Inaddition, I tried to download sequences directly from ncbi ftp site.

https://ftp.ncbi.nlm.nih.gov/refseq/release/README

wget ftp://ftp.ncbi.nih.gov/refseq/release/bacteria/
generates ----> index.html file

The html file contains all the ftp links related to the bacterial genomes.

How do I download the sequences using the index.html file ?
Is their any easy way to download the bacterial, viral genomes ?

Index of /refseq/release/bacteria on ftp.ncbi.nih.gov:21

Index of /refseq/release/bacteria on ftp.ncbi.nih.gov:21

2022 May 05 22:46 File bacteria.1.1.genomic.fna.gz (112297068 bytes)
2022 May 05 22:46 File bacteria.1.genomic.gbff.gz (90133141 bytes)
2022 May 05 22:46 File bacteria.10.1.genomic.fna.gz (2647137 bytes)
2022 May 05 22:46 File bacteria.10.genomic.gbff.gz (2242108 bytes)
2022 May 05 22:45 File bacteria.100.1.genomic.fna.gz (94215320 bytes)
2022 May 05 22:46 File bacteria.100.genomic.gbff.gz (82718913 bytes)
2022 May 05 22:46 File bacteria.1000.1.genomic.fna.gz (120745492 bytes)
2022 May 05 22:46 File bacteria.1000.genomic.gbff.gz (104026228 bytes)
2022 May 05 22:46 File bacteria.1001.1.genomic.fna.gz (120374412 bytes)
2022 May 05 22:46 File bacteria.1001.genomic.gbff.gz (102935563 bytes)
2022 May 05 22:46 File bacteria.1002.1.genomic.fna.gz (116037436 bytes)
2022 May 05 22:46 File bacteria.1002.genomic.gbff.gz (100084701 bytes)
2022 May 05 22:45 File bacteria.1003.

[1]: https://github.com/kblin/ncbi-genome-download

Answer: There is always a way directly from the FTP site:

https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/940/077/525/GCF_940077525.1_ASM94007752v1/

You likely need the file that ends in `genomic.fna.gz`.

You may want to try this program:

https://github.com/pirovc/genome_updater

Beware that going with 16 parallel downloads will likely cause NCBI to throttle down your IP number as you may be considered to abuse their resources. Sometimes going with 4-8 parallel downloads will get the job done faster than with 16.
Question: Hello ,

I have a Hiseq-4000 run data.Out of eight lane , 1,2,3,5,6,7,8 are eight index length and lane 4 is six lane index. I have made two sample sheet like sample sheet with 8 index barcode and sample sheet with 6 index. First I have demultiplexed the data (bcl2fast V2.16) with 8 index Sample sheet by using `--use-bases-mask Y*,I8n,Y*` base mask and got the out put .After wards I have tried to do demultiplex the data with 6 index length sample sheet by using `--use-bases-mask Y*,I6n,Y*` base mask. In Second case (6 base) got the following error

```
Error : One or more read masks do not match the number of cycles specified in RunInfo.xml, and a '*' was specified. There is no way to determine the number of cycles implied by the '*'. Base mask: y*'
```

Please can you tell me that how to use the base mask in such situation , same run with two type of index (6 and 8 length barcode).

Answer: Try `--use-bases-mask Y*, I6n*,Y*` for the second run with 6 bp index. Remember to specify a unique output location.

Manual for the latest bcl2fastq seems to indicate that it may be possible to do this via a single run by specifying the two bases mask this way

--use-bases-mask Y*, I8n,Y* --use-bases-mask 6:Y*, I6n*,Y*

This uses the second mask for lane 6 and the other mask for all other lanes.