cannin/biostars_qa huggingface dataset:
Question: The European Nucleotide Archive (ENA) requires certain datasets to be uploaded to the Webin upload area (`webin.ebi.ac.uk`) before they can be submitted via XML files that reference the uploaded file.

According to [this doc][1], one can use `lftp` to upload a data file. However, when I try to do this in practice, I find that several modifications have to be made to get anywhere near a login to the Webin upload area:

Attempting a log-in under the default settings causes a freeze at line
> `ls' at 0 [FEAT negotiation...]

michael_science ~ $ lftp webin.ebi.ac.uk
lftp webin.ebi.ac.uk:~> login myWebinAccount
Password:
lftp myWebinAccount@webin.ebi.ac.uk:~> ls
`ls' at 0 [Connecting...]
`ls' at 0 [FEAT negotiation...]

Attempting a log-in upon disabling FEAT negotiation and SSL encryption goes further, but still doesn't get the job done, as it freezes at line
> `ls' at 0 [Logging in...]

lftp -e 'set ftp:use-feat false; set ftp:ssl-allow false' webin.ebi.ac.uk
lftp webin.ebi.ac.uk:~> login myWebinAccount
Password:
lftp myWebinAccount@webin.ebi.ac.uk:~> ls
`ls' at 0 [Connecting...]
`ls' at 0 [Logging in...]

Has anyone successfully logged onto to the Webin upload area (`webin.ebi.ac.uk`) recently? If so, how?

**Edit 1**: Upon communication with the ENA helpdesk, it turned out to be a server problem by ENA, which seems to have been resolved for now.

[1]: http://ena-docs.readthedocs.io/en/latest/upload_01.html#general-instructions-for-uploading-files-using-a-ftp-client

Answer: I was able to connect to `webin.ebi.ac.uk` yesterday, using both `ftp webin.ebi.ac.uk` command or Transmit (FTP client for OSX), but the connection was not very stable. Today I found it not possible to reach it using the same methods.
Question: ![enter image description here][1]

Hi!

I usually work with genomics and transcriptomics data and have no experience with proteomics analysis. For the RNA-seq samples I am analyzing right now, there is proteomics analysis available, where they provide the list of significantly changing proteins. My job is to search for these changes on transcriptomics level. My question is, what does this semicolon notation mean in protemoics?

See the image attached. Does this mean, that the AURKA;AURKB complex has significantly changed? Or that this protein is either AURKA or AURKB but can't be properly identified?

[1]: /media/images/58744b56-0309-4714-8c37-23e38b1a

Answer: In datasets I've seen it means that the peptide cannot uniquely be mapped to a single protein but multimaps between the indicated proteins. Remember that in proteomics you first chop down the protein into peptides by enzymatic digestion (e.g. trypsin) so you do not measure the full protein sequence but peptides, and identical domains between related proteins can then cause multimapping.
Question: Hi Guys,

I'm pretty new to bioinformatics

I can't run the bigWigToBedGraph binary form UCSC binary list

bigWigToBedGraph: cannot execute binary file

What could be the problem?

My command was-

bigWigToBedGraph file.bw file.bedgraph

Is there an issue with my usage?

Answer: Not to worry - got it sorted!

Rookie mistake - grabbed the wrong file!
Question: Hi folks.

I'm fine-mapping an association to a gene using next-generation sequencing data in African-Americans. I have a list of a few thousand variants on chromosome 12 listed by position based on UCSC hg19 build 36.

I want to retrieve the reference/variant alleles and minor allele frequency from 1000 genomes project for YRI samples for comparison to my own sequencing data. How might I best do this without downloading the 1000 genomes data and recomputing allele frequencies? Is there a way to query ensembl or UCSC for this information?

Thanks in advance.

Answer: Thanks for all the responses. I actually got a very helpful response from a friend of mine. It was dead simple to download and compile vcftools and tabix on my virtual linux system. It took about 30 seconds to download my genotype calls of interest from 1000 genomes using tabix:

tabix -f -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/supporting/AFR.2of4intersection_allele_freq.20100804.genotypes.vcf.gz 12:101266000-101422000 > afr.vcf

Then less than a second to analyze this data for allele frequency using vcftools:

vcftools --vcf afr.vcf --freq --out freq-afr

Which gives me exactly what I want.

CHROM POS N_ALLELES N_CHR {ALLELE:FREQ}
12 101266049 2 348 A:0.87069 G:0.12931
12 101266292 2 348 A:0.954023 T:0.045977
12 101266307 1 10 G:1
12 101266725 1 10 T:1
12 101266729 2 348 G:0.991379 A:0.00862069
12 101266756 2 348 T:0.985632 G:0.0143678

Now, I have frequency info for my data mapped to UCSC hg19/b36. Now I just need to figure out how to match these up to the 1000 genomes frequencies using their positions.
Question: Hi,

I've used freebayes for calling variants on a set of bamfiles with 1 sample in each bam file.
I want to somehow get a good overview over the variants that also non-bioinformaticians can look at.

The VCF file looks something like this:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1 SAMPLE2 SAMPLE3 SAMPLE4 CONTROL
chr1 46 . ACACC CCACA 8.58901e-15 . AB=0;ABP=0;AC=0;AF=0;AN=5;AO=2;CIGAR=1X3M1X;DP=137;DPB=147;DPRA=0.314961;EPP=3.0103;EPPR=17.5948;GTI=0;LEN=5;MEANALT=1;MQM=33;MQMR=54.4254;NS=5;NUMALT=1;ODDS=78.0575;PAIRED=1;PAIREDR=0.977612;PAO=0;PQA=0;PQR=506;PRO=17;QA=14;QR=4464;RO=134;RPP=7.35324;RPPR=77.9423;RUN=1;SAF=1;SAP=3.0103;SAR=1;SRF=83;SRP=19.6042;SRR=51;TYPE=complex;technology.ILLUMINA_HiSeq2500=1 GT:DP:RO:QR:AO:QA:GL 0:81:80:2614:0:0:0,-10 0:6:6:218:0:0:0,-10 0:10:8:262:2:14:0,-10 0:19:19:671:0:0:0,-10 0:21:21:699:0:0:0,-10
chr1 2146 . G A 9.16135e-15 . AB=0;ABP=0;AC=0;AF=0;AN=5;AO=28;CIGAR=1X;DP=251;DPB=251;DPRA=0;EPP=5.80219;EPPR=11.1996;GTI=0;LEN=1;MEANALT=1;MQM=18.9643;MQMR=57.3318;NS=5;NUMALT=1;ODDS=92.9516;PAIRED=0.964286;PAIREDR=0.977578;PAO=0;PQA=0;PQR=0;PRO=0;QA=633;QR=8052;RO=223;RPP=63.8115;RPPR=5.82445;RUN=1;SAF=11;SAP=5.80219;SAR=17;SRF=99;SRP=9.09627;SRR=124;TYPE=snp;technology.ILLUMINA_HiSeq2500=1 GT:DP:RO:QR:AO:QA:GL 0:145:136:4920:9:154:0,-10 0:27:22:788:5:129:0,-10 0:16:13:471:3:75:0,-10 0:13:9:333:4:100:0,-10 0:50:43:1540:7:175:0,-10

I've tried using both gatk and snpsifts variantotable/extractfields tools, but it will only sum all the INFO counts and not keep them separated.

My goal is to get something like this:

CHR POS EFF.FIELD1 REF ALT QUAL SAMPLE1.RO SAMPLE1.AO SAMPLE2.RO SAMPLE2.AO SAMPLE3.RO SAMPLE3.AO SAMPLE4.RO SAMPLE4.AO CONTROL.RO CONTROL.AO
chr1 46 "geneX" ACACC CCACA 8.58901e-15 80 0 6 0 8 2 19 0 21 0

But the current options would give something like:

CHR POS EFF.FIELD1 REF ALT QUAL RO AO
chr1 46 "geneX" ACACC CCACA 8.58901e- 134 2

My feeling is that this should be fairly common, so just want to make sure I'm not reinventing the wheel here.
Anyone have any idea?

Answer: You should use the latest version of bcftools, and you should get exactly what you want. However you would have to take time to make the `-f` string perfect. For example:

bcftools query --print-header -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT\t%RO\t%AO]\n' file.vcf.gz
Question: Hello everyone!

I am trying to learn more about the SAM/BAM format and how reads are mapped to the genome.

I understand that reads can be mapped to the forward (sense) and reverse (antisense) strand.

I understand that in every read pair there is a first-in-pair and a second-in-pair read - and that this probably has something to do with the order in which they came off the sequencer and not the position in the genome? (but maybe im wrong)

I understand that weird things can happen. A pair of reads can map millions of base pairs apart. Both reads can map to the same strand. But can reads map in both upstream and downstream direction - i.e. the read sequence once reversed (but not reversed complimented) was mapped.

With three variables thats 32 different combinations for paired reads (and another 4 for singletons).

With two variables, there's only 8 different combinations (and 2 for singletons)

Regardless, does anyone know of any tool available which will give you the statistics for the 32+4 or 8+2 combinations in your BAM file? If not, i will add it to the new version of metaflagstat.

Here's a picture of what i'm trying to explain :)

Answer: Using proper nomenclature, a given read can map to the + strand or reverse complemented to the + strand...which is the same as mapping to the - strand. It will never be mapped reverse to anything. Things a bit more complicated if you're talking about bisulfite converted reads, but aside from that those are the two possibilities. This is at least true with current technologies, but will likely remain so in the future due to the energetics of 3'->5' base addition.

With PE reads, any relative orientation of the mates to each other is possible, though some are vastly more likely than the others.
Question: How to know a RNA-seq is polyA or total-RNA-seq in GEO?

Does GEO require that kind of meta data during submission?

Thanks,

Answer: As a rule of thumb, if a project doesn't explicitly mention looking at total RNA then the odds are high that it's polyA selected. When in doubt, you can always align against rRNA, since the alignment rate for polyA selected samples will be very low (<<1%), while rRNA-depleted (aka "total RNA") will have >1% rRNA (sometimes up to ~30%).
Question: Hi there,

I always used TCGAbiolinks to get raw count for TCGA projects like below:

expquery <- GDCquery(project = "TCGA-KIRC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
GDCdownload(expquery,directory = "GDCdata")
expquery2 <- GDCprepare(expquery,directory = "GDCdata",summarizedExperiment = T)
expMatrix <- TCGAanalyze_Preprocessing(expquery2)

However, it does not work today and it seems there is no HTSeq - Counts

Error in GDCquery(project = "TCGA-KIRC", data.category = "Transcriptome Profiling", :
Please set a valid workflow.type argument from the list below:
=> STAR - Counts

Therefore, I used STAR - Counts to download the data, which has completely different format of what I downloaded before using HTSeq - Count. The expression matrix for each sample has more columns including fpkm_unstranded and tpm_unstranded.

Actually the data downloaded using STAR - Counts is much more useful but I do not know how to extract the files to a readable expression matrix (ENSEMBL ID as rownames and TCGA tumor barcode as column names) because the GDCprepare() function fails to work on it:

expquery2 <- GDCprepare(expquery,directory = "GDCdata",summarizedExperiment = T)
| | 0%
Error in readr::read_tsv(file = f, col_names = TRUE, progress = FALSE, :
unused argument (show_col_types = FALSE)
Error in if (value == n) { : argument is of length zero

Anyone has any solutions? Many thanks in advance!

Answer: I have found the solution!
Just update the TCGAbiolink package from the Github project, here is the link:
https://github.com/BioinformaticsFMRP/TCGAbiolinks/issues/495

I had some problems during the installation, I had to reinstall some packages to make it work, here they are:

library(GenomicDataCommons)
library(GenomeInfoDbData)
library(SummarizedExperiment)
library(ExperimentHub)
library(AnnotationHub)
Question: Im using this code to make based on log2foldchange and padj value ,im getting the plot but i want those value for my reference how do i extract the same .

alpha <- 0.05 # Threshold on the adjusted p-value
cols <- densCols(res$log2FoldChange, -log10(res$pvalue))
plot(res$log2FoldChange, -log10(res$padj), col=cols, panel.first=grid(),
main="Volcano plot", xlab="Effect size: log2(fold-change)", ylab="-log10(adjusted p-value)",
pch=20, cex=0.6)
abline(v=0)
abline(v=c(-1,1), col="brown")
abline(h=-log10(alpha), col="brown")

gn.selected <- abs(res$log2FoldChange) > 2.5 & res$padj < alpha
text(res$log2FoldChange[gn.selected],
-log10(res$padj)[gn.selected],
lab=rownames(res)[gn.selected ], cex=0.4)

when i view gn.selected i get only logical value that is true or false

Any help or suggestion would be highly appreciated

Update I'm doing this

> DF <- DF[DF$log2FoldChange > 1.5 & DF$padj < 0.05,]

is that suffice and am i doing it correctly ?

Answer: Easy Volcano plot for gene expression data in Python
https://reneshbedre.github.io/blog/volcano.html
Question: Hi all,

I've designed primers to identify a strain of interest.

Just to be sure that the primers matche only the genome of this strain and no other genomes, I am using blastn analysis.

My input file is :

>Seq_1
ATGTGCTAGCTATGGCG
>Seq_2
GTCGGTAATCGATTCGA

And my database is a multifasta file containing several genomes including the genome of my strain of interest.

My command line is the following :

> blastn -query Sequences.fasta -db genomes.fasta -out Results_blastn.txt -outfmt 6 -evalue 0.1

**The results file is blank.**

I do not understand how I've got no results even on the genome of my strain of interest ...

Do someone have an idea on how I should proceed to fix this issue?

Thank you so much in advance for your much appreciated help!

Audrey

Answer: Just add `-task blastn-short` to your command. By default, blastn uses megablast (optimized for highly similar long sequences).
Question: When working with alternative transcripts, I see a decimal ending to the transcript id. What does the decimal value add to the transcript meaning? For example, a mutation, BRAF p.V600M is annotated by VEP with:

csn="ENST00000479537.1:c.82G>A_p.Val28Met"

hgvsc="ENST00000479537.1:c.82G>A"

transcript_id="ENST00000479537"

&

csn="ENST00000288602.6:c.1798G>A_p.Val600Met"

hgvsc="ENST00000288602.6:c.1798G>A"

transcript_id="ENST00000288602"

I see the transcript IDs for BRAF listed at Ensembl as they are given in the "transcript_id" fields above. What does the ".1" and ".6" mean for each transcript above?

Answer: They are version numbers. They're not useful for most practical purposes because if the sequence changes significantly, then the ID itself is changed. See here for an explanation of versioning in Ensembl.

EDIT: update link
Question: Dear Biostars, Hi

I have a list of **zebrafish** TFs (both family and members) upregulated in my condition1 and same data for condition2. (I also have a list of all transcription factores showing blastx hit with all my RNA-seq analysis)
I want to show them in a network (or maybe two networks).

I saw a [similar network][1] created using *Cytoscape*.

Could you please kindly help me in this issue? how to ready files as *cytoscape* input? how to create circular network showing each DE TFs in different colour?

or even another package/website for doing this ?

Thanks

NOTE: e.g. of TFs families Up-regulated in my condition 1:

BX649453.2

baz2a

gtf2ird1

tcf3b

mef2d

lhx2b

PBRM1

cers5

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4428447/figure/F6/

Answer: Hey Farbod,

You first need to decide whether to make a protein interaction network or co-expression network.

Then you can use [STRING][1], [GeneMANIA][2] or [BioGRID][3] to build your network.

Once the network is built, you can import it into Cytoscape and change the layout to circle it, and look the same as the picture you showed.

To change the color or shape of only a few nodes, you can add a new column to your network in cystoscape and for example enumerate only the node of interest, there you can only change the shape and color of that node without changing the other.

Best Regards,
Leite

[1]: https://string-db.org/
[2]: https://genemania.org/
[3]: https://thebiogrid.org/
Question: For one of my samples, I have two runs.

sample1
+-- run1/
| +-- *.bax.h5
+-- run2/
+-- *.bax.h5

My workflow involves using **bax2bam** to convert *.bax.h5* to *unaligned bam*. Then using **pbalign** to convert *unaligned bam* to *aligned bam*.

How do I combine these two runs into one sample? And at what point in my workflow do i do that?

Answer: I managed to get it to work. Not sure if this is the best way.

I converted `*.bax.h5` to BAM separately for each run.

bax2bam -o "sample1_run1" $(find . -name "*.bax.h5")
bax2bam -o "sample1_run2" $(find . -name "*.bax.h5")

This creates scraps and subreads BAM files. I merged the subreads BAMs.

samtools merge sample1.bam sample1_run1.subreads.bam sample1_run2.subreads.bam

Then created the **.pbi** index file using SMRT/5.0.1

pbindex sample1.bam
Question: The sequencing of this ChIP-seq H3K4me3 experiment (human, single-end 50bp) has a problem with the number of reads per sample. The person that did the sequencing did not arrange correctly the samples in the lanes, causing the bias you see in the picture (first, upper one). I am trying to "rescue" it doing downsampling, but wanted to make sure I am doing it correctly without loosing too much coverage that I loose robustness in whatever finding I have later on.

I did the downsampling as described in Ming Tang's blog : https://crazyhottommy.blogspot.com/2016/05/downsampling-for-bam-files-to-certain.html . Here there is an example line of code:

samtools view -s 3.6 -b my.bam -o subsample.bam

I did three downsamplings to see what happens: one to 20 M reads, one to 25M and one to 10M. In the figure below you only see raw data, 25M and 10M in that order (from top to bottom). My questions are:

1) How much you think I can downsample without having an unreasonable coverage? Is 10M ok or given that is single-end more should be used? (for example, in Girdhar et al 2018's schizophrenia chip-seq paper in humans they used 13M but it was paired-end sequencing).

2) If I end up deciding 10M is too little and I go with the 25 or 20M one (please, look at the second plot), you can still see a difference between the samples. Is that difference small enough to continue analysing or you think that's already too much?

3) In case in 2) you consider this difference is too big and that that's not good, you have another suggestion of what I could do?

Thank!

![enter image description here][1]

[1]: https://i.ibb.co/dbfCgTs/Screenshot-2019-06-28-at-11-01-42.png

Answer: > 1) How much you think I can downsample without having an unreasonable
> coverage? Is 10M ok or given that is single-end more should be used?
> (for example, in Girdhar et al 2018's schizophrenia chip-seq paper in
> humans they used 13M but it was paired-end sequencing).

The current ENCODE standard for H3K4me3 is a minimum of 20M mapped and filtered reads, however the previous standard was 10M. Because the important aspect of ChIP-seq analysis is the location of the alignments, not the actual sequence content of the reads, as long as the alignment rate is good with single-end reads then it's not much different from a similar paired-end depth.

> 2) If I end up deciding 10M is too little and I go with the 25 or 20M
> one (please, look at the second plot), you can still see a difference
> between the samples. Is that difference small enough to continue
> analysing or you think that's already too much?

It depends on the analysis you will be doing with the data. If your analysis is going to be based on looking at peak calls, for example, I would go ahead and do your peak calling on the full data and see how affected that is by the differences in depths.

Btw, are the numbers you're showing the raw reads or the number of reads after mapping and removing duplicates?
Question: Is there anyway to get BioGrid gene interaction network automatically using R?

Answer: Coincidence! I've just been putting together an R client for BioGRID. Very rough initial version is at https://github.com/rajarshi/cdkr/tree/master/rbiogrid. To install

```
library(devtools)
install_github("cdkr", "rajarshi", subdir="rbiogrid")
```

You will need to get an access key from the BioGRID website. Right now it has no tests and has basic functionality to get interactions by gene names. It sets a number of arguments by default (e.g., `selfInteractionsExcluded = TRUE, includeInteractorInteractions = FALSE`) - you can view the source of the function to see them.

Example usage is

```
library(rbiogrid)
set.access.key('XXXX')
get.interactions('PLK1')
```
Question: I've assembled a known E. coli genome only using nanopore reads. I want to check if the genes predicted matches the ground truth. I tried Prokka but it predicts way too many CDS (~10000). Although most of them are 'hypothetical proteins' so I'm guessing those are false positive ORFs caused by frameshifts as a result of nanopore's systemic indel errors.

Are there programs that work well specifically with nanopore erros in mind? Glimmer and GeneMark both have long histories but I'm guessing they're optimized for short reads.

Answer: Basically what you got is a very bad assembly, and there's no reason to look at gene predictions. There are two possible things you can do here

- if you want a good assembly, you also need to polish yours with Illumina reads (or do a hybrid assembly with e.g. Unicycler) - that would ensure you have the least errors;
- if you just want to get the idea which genes are present, take a reliable protein database (e.g. 18k well curated proteins used by Prokka) and run tblastn on your genome assembly. You won't be able to distinguish genes/pseudogenes, but otherwise you'll get a pretty good idea about gene presence/absence.
Question: I have a multi-fasta file containing sequences with headers such as:

```
>AAST01014508.1|(1..6240)|LTR/Pao|ROOA_I-int:ROOA_LTR
>AAST01026747.1|c(2745..6820)|LTR/Pao|ROO_I-int
```

In the second section of the header where the position is given I want to index all of the sequences using this and print to file all those sequence from the C strand, e.g. `>AAST01026747.1|c(2745..6820)|LTR/Pao|ROO_I-int`

I was able to create a function easily which did this:

```
def get_compstrandTE(record):
parts=record.id.split("|")
assert len(parts) ==4
return parts[1]
```

However, I am now stuck as to how to search through the keys in the dictionary and find only those containing the 'c' and write those to a file. I tried using the example from the BioPython manual but kept running into difficulties.

If anyone has any suggestions I would really appreciate it.

Answer: You don't need a dictionary for this task.

Something like this should work, loosely based on of the filtering examples from [the Biopython Tutorial][1] using a generator expression:

```
from Bio import SeqIO
input_file = "big_file.fasta"
output_file = "complements.fasta"

def wanted(record):
"""Returns True if name scheme suggests from complement stand."""
parts = record.id.split("|")
assert len(parts) == 4
return parts[1].startswith("c")

records = (r for r in SeqIO.parse(input_file, "fasta") if wanted(r))
count = SeqIO.write(records, output_file, "fasta")
print("Saved %i records from %s to %s" % (count, input_file, output_file))
```

[1]: http://biopython.org/DIST/docs/tutorial/Tutorial.html