cannin/biostars_qa huggingface dataset:
Question: Question is in the title.
Imagine you have many bigwig track...and you want to display the same scale for each tracks but you don't want to change viewlimits one by one for all of them.
Is there a way to do that in a one time shot for all using maybe session file or Supertrack ,multiwig, composite arguments.

Answer: You can check the "Oct. 11, 2019 New "group auto-scale" option" [here][1].

Example track file for group of tracks:
[1]: http://genome.ucsc.edu/goldenPath/newsarch.html#101119

bigWig

track myComposite
compositeTrack on
shortLabel test
longLabel my_comp_testing
type bigWig

track myFirstSubtrack
parent myComposite on
shortLabel test_1_short
longLabel tes_1_long
type bigWig
visibility full
maxHeightPixels 70:70:32
configurable on
autoScale off
viewLimits 0.0:1.0
alwaysZero on
priority 0.1
bigDataUrl file1.bw
color 0,0,0

track mySecondSubtrack
parent myComposite on
shortLabel test_2_short
longLabel tes_2_long
type bigWig
visibility full
maxHeightPixels 70:70:32
configurable on
autoScale off
viewLimits 0.0:1.0
alwaysZero on
priority 0.1
bigDataUrl file2.bw
color 0,0,255
```
Question: So I'm trying to integrate two datasets together for 10 different samples (taken in different batches). To mitigate for that I looked into Harmony but was not sure how to essentially apply Harmony from two datasets and get fragment files for later analysis in ArchR.

Any suggestions or pointers?

Answer: `ArchR` has built-in support for Harmony correction;

```
projHeme2 <- addHarmony(
ArchRProj = projHeme2,
reducedDims = "IterativeLSI",
name = "Harmony",
groupBy = "dataset"
)
```

https://www.archrproject.com/bookdown/batch-effect-correction-wtih-harmony.html
Question: Two map/ped files refused to merge for some reason.

I followed this method here (https://www.biostars.org/p/56640/ by [ff.cc.cc][1]) to reduce the datasets to the same list of sites.

I used

plink --bfile Deven --bmerge Harvard.fam Harvard.bim Harvard.bed --merge-mode 6 --out mismatchB

to identify the sites that had flipped genotypes, and I straightened them out using

plink --bfile Harvard --flip mismatchB.missnp --recode --make-bed --out HarvardFlipped

I converted them to map/ped from the binary format, and I tried to merge them and got an inexplicable map file error

```
ERROR: Problem with MAP file line:
HGDP00001 HGDP00001 0 0 1 -9 G T A A C C G G A G G G A A C T G G A G T T T T G G C C C C A A T C C C C C G G G G T T T T G G A G T T C C A A A G C T A G C C G G A A G G C C T C C T A G A G G G G G A A T C T T T T G G G G A G G G C C T C G A G G G G C C T
```

I looked at the map and ped files and everything was fine. Any clue what is going on? I know that the Harvard dataset originally has coordinates using Hg18 and the Deven dataset has coordinates the using Hg19 (but the rs ids come from the same annotation file, which provided the conversion of Affy #s to dbSNP #s)

[1]: https://www.biostars.org/u/4466/

Answer: UPDATE: I merged them as binary, and then converted them to map/ped and it worked... You all can probably ignore this.
Question: Hi,

Easy question here. I am split up my vcf file up by chromosome to try and save some time. I want to combine the files now to have one complete vcf file with all the chromosomes. Do I wanted to merge the files or do I want to concatenate the files?

I just ran the command

cat *.vcf > cat.chr1.chr2.vcf

This was only on two of the chromosomes and when I went to look at the files, there was only data for the first chromosome and nothing for the second. Am I doing something wrong here? I could also try vcf tools if this "cat" command wont work for what I am trying to do.

Answer: first alternative, use one vcf file to get the header, and the concatenante all the other vcf without the header (this was your error)

grep '^#' chr1.vcf > merge.vcf
grep -v '^#' chr1.vcf chr2.vcf chr3.vcf chr4.vcf >> merge.vcf

second alternative: use picard gatherVcfs http://broadinstitute.github.io/picard/command-line-overview.html, which is going to check the headers, orders etc...

java -jar picard.jar GatherVcfs I=chr1.vcf I=chr2.vcf I=chr3.vcf I=chr4.vcf O=merged.vcf
Question: Good morning,

I want to analyse public SHort Read archive from an RNA-seq experiment.
I am trying to quantify them using sailfish. it is needed to provide the library type as an argument : -l ""
But I only that my data are paired-end ( dataset GSE89063)
How can I know what to write for this argument ?
Thanks a lot

Have a good day

Answer: Hi,
You can also check by performing some steps to know the strand specificity, please check the link below.

http://onetipperday.sterding.com/2012/07/how-to-tell-which-library-type-to-use.html

Hope it works with you.
Question: Dear all,

I have FASTQ files and on start of my read I have 7 nucleotides tag - I would like to extract reads with this specific tag.

I would like to search in first 15 nucleotides of my reads, if match - extract this read to new fastq files.

Thank you for any ideas or help.

Answer: In order to stay old-school, you can use the FastX toolkit's [barcode splitter][1].

You need a text file with your tag (let it be myTag.txt) and then you can run (with N mismatches):

cat sequence.fastq | /usr/local/bin/fastx_barcode_splitter.pl -Q33 --bcfile myTag.txt --bol --mismatches $N --prefix out --suffix .txt

In order to remove the tag-sequence, you can use the fastx_trimmer from the same toolkit

Using the Fastx-tools, do not forget to add the -Q33 option.

[1]: http://hannonlab.cshl.edu/fastx_toolkit/commandline.html#fastx_barcode_splitter_usage
Question: I want to store a docker image with the input data, code, software, intermediate results and results. This will be a massive file.

Is there something like the sequence archive for people who want to upload a docker image of their whole study? If not, what is the best/cheapest way to store such a file?

Any related thoughts not directly answering the question are also welcome.

Answer: Here are a few options that come to my mind:

- Host the data on a public repository, for example [synapse][1], which has been created for this exact purpose. You can keep the data private during the publication, and then release it afterwards. Hosting is free of charge.

- As previously noted, the docker image should be small. In the docker image, download the data from synapse or wherever you host the data, and pull the code from github or bitbucket

- if you really need to host large datasets, bitbucket has no space limits (although large repositories may incur into limits)

- I remember some bioinformatics-related journals advertising that they didn't have limits for the size of supplementary materials posted. This is to encourage people to upload all the relevant data and facilitate reproducibility. However right now I don't remember the name of the journals using this policy. Sorry :-)

[1]: https://www.synapse.org/
Question: Hello everyone!

*Today I was trying to check whether all SNPs in my list (thousands of SNPs, GRCh37) belong to the particular chromosome (let's say chr 19) but I had some issues with it. I am doing everything on R and my strategy was this:*

1. To prepare list of SNP as a file (all starts with rs...)
2. To download the human chromosome 19 from dbSNP
3. To search my list as a pattern within downloaded chromosome

However, downloading process of chr 19 is going on for many hours already... And the fastest download option doesn't give the proper information for analysis.

I would love to know whether there is another way to check the list of SNPs for belonging to the particular chromosome in R?

Thanks in advance!

Answer: Usign *biomaRt*:

#set up
library(biomaRt)
snpmart = useEnsembl(biomart = "snp", dataset="hsapiens_snp")
#Other useful functions
# listAttributes(snpmart)
# listFilters(snpmart)

Get SNPs that match with our list of IDs, here we have 2 SNP from chr1 and 1 SNPs from chr2:

getBM(attributes = c("refsnp_id", "chr_name", "chrom_start", "chrom_end"),
filters = c("snp_filter"),
values = list(c("rs17599629", "rs1218582", "rs11902236")),
mart = snpmart)

# refsnp_id chr_name chrom_start chrom_end
# 1 rs11902236 2 9977740 9977740
# 2 rs1218582 1 154861707 154861707
# 3 rs17599629 1 150685811 150685811

Now, to return only the ones which are on chr2, we can add `"chr_name"` filter:

getBM(attributes = c("refsnp_id", "chr_name", "chrom_start", "chrom_end"),
filters = c("chr_name", "snp_filter"),
values = list(c(2),
c("rs17599629", "rs1218582", "rs11902236")),
mart = snpmart)

# refsnp_id chr_name chrom_start chrom_end
# 1 rs11902236 2 9977740 9977740
Question: Say I have two VCFs with 100 samples in each file. Each VCF was joint-called separately and now I want to merge the variant calls.

Do I need to left-align the INDELs in the merged VCF? I've used `bcftools norm` in the past and got odd results. It seems that `vt` is a better tool for this.

Is left-aligning only useful for common variants? If I'm interested in rare variants (<0.5% AF) would left-alignment actually matter?

Thanks

Here's an example of bcftools norm

Original VCF

chr7 157009949 . AGCGGCGGCGGCG AGCGGCGGCGGCGGCGGCG,A,AGCGGCGGCGGCGGCGGCGGCG,AGCGGCGGCG,AGCGGCGGCGGCGGCGGCGGCGGCG,AGCGGCGGCGGCGGCG

---

Left-Aligned VCF (with multiallelics split into biallelic calls)

chr7 157009949 . A AGCGGCGGCG
chr7 157009949 . A AGCGGCGGCGGCG
chr7 157009949 . A AGCGGCG
chr7 157009949 . A AGCG
chr7 157009949 . A AGCGGCGGCGGCGGCG
chr7 157009949 . AGCG A
chr7 157009949 . AGCGGCG A
chr7 157009949 . AGCGGCGGCGGCG A

---

Answer: Suppose one of your VCF files has a non-left-aligned insertion, represented as REF=AG, ALT=AGT, starting at position 99999, and another file has an insertion represented as REF=G, ALT=GT, starting at position 100000. If you don't left-align, these may not be recognized as the same variant, and downstream analysis will suffer.
Question: Dear all,

for a list of Structural Variants (including deletions, duplications, inversions, translocations), either in VCF or BEDPE format, we would like to have the gene annotations, and the lists of the following sets of genes :

-- fusions (if both breakpoints are in exons, introns, utrs)
-- truncations (if only one breakpoint is in exon, intron, utr; and the other breakpoint is in intergenic area)
-- the genes in the areas that are deleted, duplicated, inverted

Although I wrote some scripts in perl based on Annovar , thought that we could get all these annotations with a package that is already available ?

thanks a lot,

-- bogdan

Answer: SVs are problematic for many pipelines/software as, unlike SNVs and small indels, each event involves at least two genomic loci.

Be aware that not all callers correctly classify events. Many callers will classify events purely on their break-end position and orientation. This results in deletion calls even when there is no copy number change to support the event (most callers), or an inversion calls even when only one of the two inversion breakpoints actually exist (e.g. DELLY). For simple germline analysis this is probably ok, and you can just ignore all large or inter-chromosomal events but for highly rearranged genomes (eg cancer), things are much more complicated.

> thought that we could get all these annotations with a package that is already available

What you're asking is really two separate processes: one for looking at the intervening sequence of simple events, and another for break-end overlap for fusions/interchromosomal/complex events.

If you're familiar with BioConductor then you can do the first part relatively easily for a BEDPE: just convert to GRanges intervals and calculate overlaps against the BioConductor annotation package for your organism.

For the second part you might be interested in my [StructuralVariantAnnotation][1] package. It's key feature is conversion of VCFs generated by a number of popular SV callers into a GRanges object containing break-end coordinates. Once in GRanges format, you can again use the BioConductor annotation packages to calculate feature overlap.

[1]: https://github.com/PapenfussLab/StructuralVariantAnnotation
Question: Here's another beginner BioPython question from me...

I'm running some genome assemblies for someone who has some new Illumina sequence data and also had done some sequencing a few years ago. They have some Sanger and 454 sequences (a couple thousand sequences with a couple thousand base pairs for each) and everything is in fasta format. They are trusted contigs and I'd like to use the data to help bridge assembly gaps of Illumina sequencing data (from the same extracted DNA that was kept frozen for a few years -- same bacterial strain). There are no quality files for the older data and if they ever existed they are long gone from the data directories.

What I want to do is create FASTQ files for the old Sanger and 454 FASTA files with a phred score of 40 -- I think that would be a `I` under PHRED-33.

(I know there are assemblers out there that take fasta files inputs along with fastq -- but my in-house makefile and pipeline does not and I really would like to incorporate these sequences in our current workflow without a lot of changes).

There are lots of resources for converting from fasta to fastq *WITH* a qual file -- see [here](/p/16458/), [here](/p/85929/), [here](http://qiime.org/scripts/convert_fastaqual_fastq.html), [here](http://seqanswers.com/forums/showthread.php?t=16925), and [here](http://sequenceconversion.bugaco.com/converter/biology/sequences/fasta_to_fastq-solexa.php); but I've looked all over but can't find anything for converting without a qual file.

[Here is a public perl script](https://code.google.com/p/fasta-to-fastq/) that does exactly what I want it to do -- but I can't get it to work for me and I keep getting an error.

As an exercise, I quickly wrote the script below, but I can't get it to provide the quality score as the length of the sequence. The quality string needs to be the same length as the input sequence string. Anyone have any hints for me? Any help is appreciated.

Here's my code so far:

#!/usr/bin/env python

"""
Convert FASTA to FASTQ file with a static quality read

Usage:
$ ./fasta_to_fastq NAME.fasta NAME.fastq
"""

# import libraries
import sys, os
from Bio import SeqIO
from Bio.SeqIO.QualityIO import PairedFastaQualIterator

# Get inputs
fa_path = sys.argv[1]
fq_path = sys.argv[2]
phred_quality = "I"

# Check inputs
if not os.path.exists(fa_path): raise Exception("No file at %s." % fa_path)

# make fastq
with open(fq_path, "w") as handle:
records = PairedFastaQualIterator(open(fa_path), print("phred_quality"))
count = SeqIO.write(records, handle, "fastq")

# print fastq
print "Converted %i records" % count

Answer: Since Matt & Istvan have given you SeqRecord based solutions using SeqIO, here's a light weight solution using strings:

from Bio.SeqIO.FastaIO import SimpleFastaParser
with open("input.fasta") as in_handle:
with open("output.fastq", "w") as out_handle:
for title, seq in SimpleFastaParser(in_handle):
out_handle.write("@%s\n%s\n+\n%s\n" \
% (title, seq, "I" * len(seq)))

Note that faking quality scores like this is probably not a good idea, and to be used with caution.
Question: Hi, I am completely a newbie in the field of NGS-RNAseq. But I have a set pipeline to do RNAseq analysis. But to start with my analysis I need to go with TCGA database. But I cannot find any fastq files and the normal samples for the respective data set. Please guide me how to work with TCGA database to download fastq files for specific type of cancer with the respective normal samples to do the RNA-seq analysis?

I've queried [Breast invasive carcinoma][1] [BRCA], which gave me 1094 mRNA cases. And number of NT and TN samples are listed out. but for RNAseq level 1 no data is available, for level 2 only vcf files are available and for level 3 all the outputs and analysis files are available.

Thanks in advance

[1]: https://tcga-data.nci.nih.gov/tcga/tcgaCancerDetails.jsp?diseaseType=BRCA&diseaseName=Breast%20invasive%20carcinoma

Answer: The Level 1 data are controlled access for ethical reasons. See that "Controlled Access Data" section on this site for details and the process.

https://tcga-data.nci.nih.gov/tcga/tcgaAccessTiers.jsp
Question: So, I have a list of variants I want to annotate based off several databases in .bed format. The rest of the databases I want to include are off of ENCODE, but despite my searches I cannot for the love of my life download the data off ENCODE based off cell line.

I always end up at [this page](http://genome.ucsc.edu/ENCODE/downloads.html), which provides me ALL the databases, but not by cell line. I found a search function [here](http://www.epigenomebrowser.org/cgi-bin/hgFileSearch?), but for some reason, whenever I run a search, it doesn't return any results.

I am looking for .bed files listing sites of DNA methylation, DNase-seq, Chip seq (TFBS and histone modifications) for HepG2, K562, GM12878, and Human Skeletal Muscle cell lines each.

I am basically asking: How can I get these off ENCODE? Thanks.

Answer: You go to the Encode Data Matrix for hg19 -> Click on the green boxes -> You will see a lot of tracks under the column name Track Name -> There is a small arrow at the end of the tracks which will expand the track, click on it -> You will get the links to `narrowpeaks.bed.gz` or `broadpeaks.bed.gz` which you can directly download. You might not get every data that you need though.
Question: I'm studying a species where there is a reference genome that is assembled *only* at the scaffold level ("unplaced scaffolds"). See here https://www.ncbi.nlm.nih.gov/genome/annotation_euk/Geospiza_fortis/101/.

My question is

- Do people generally treat a reference genome at the scaffold level as
if each scaffold would be a chromosome?
- Should a scaffold vs
chromosome level reference genome be treated differently?
- What are
the main challenges for using a reference genome that is only at the
scaffold level?

Basically, I often read in population genetics textbook that we have to study "chromosome". But I have a hard imagine, when having only scaffold, how the theory applies differently.

Answer: All bioinformatic applications, including mapping (to a FASTA), will proceed as usual. For annotation, you'll also be fine as long as the contig/scaffold names are the same as the annotation file you're using. This isn't a major problem - the human reference, for example, has a number of unlocalized scaffolds and patch scaffolds that are relevant for annotation but aren't in the set of more thoroughly characterized chromosomes. Contigs/scaffolds/chromosomes are often treated the same by most applications.

So, to answer your questions:

1. More-or-less.
2. Not really, but realize that your analysis may be impacted if all the scaffolds can't be localized to true chromosomes (e.g., are two scaffolds in LD because they're next to each other in reality?).
3. Outlined above, but they're generally captured under issues resulting from 'assembly uncertainty' sensu lato.
Question: Hi,

I would like to analyse the content of a GTF file. I am quite able with R and dplyr, so I would like to transform my GTF file into a data frame to facilitate my analysis. Does anybody know of any tool to do this?

Thanks.
Best, C.

Answer: try `gtf_df=as.data.frame(gtf)` after importing via import function from rtracklayer.

The code would be:

```r
gtf <- rtracklayer::import('celegans.gtf')
gtf_df=as.data.frame(gtf)
```
Question: I am new to gene ontology and was wondering the best tool to use to test if a gene set is enriched for certain GO terms in a non-model prokaryote?

Answer: Not a lot is out there regarding gene set enrichment in Prokaryotes.

A few links that might be helpful and give you a push:

- https://www.biostars.org/p/85909/
- https://www.biostars.org/p/50960/

An excerpt from **Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists**

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/

> Several enrichment tools were specifically designed for these less popular species, such as WEGO for rice ([54][1]); easyGO for crops ([66][2]); FINA for prokaryotes ([58][3]); CLENCH for Arabidopsis ([21][4]); JProGo for prokaryotes ([48][5]); BayGo for *Xylella fastidiosa* ([52][6]).

[Protein Localization Analysis of Essential Genes in Prokaryotes][7]

HTH

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B54
[2]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B66
[3]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B58
[4]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B21
[5]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B48
[6]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2615629/#B52
[7]: http://www.nature.com/srep/2014/140808/srep06001/full/srep06001.html
Question: Hi,

I found [this useful website](http://cbio.mskcc.org/tcga-generanker/index.jsp) based on mutliple cancer gene list. I searched, but I cannot find any information on how this tcga-generanker is made ? and which statistics it used to compute the score ? And also it says Gene Ranker: TCGA GBM 6000. Is it only for GBM (Glioblastoma Multiforme) .

Answer: and which statistics it used to compute the score ?

If you look at their TCGA generator sources page, you can see that their list of candidate genes comes from a collection of lists of from other sources. Each list has a specific weight, and from the way it looks, it seems that they attributed scores based on the confidence they had on each list. Then, each gene is ranked with the sum of of the scores of all the lists to which it belongs

Is it only for GBM (Glioblastoma Multiforme) .?

The database contains genes specific to other tissues (e.g. COL11A1 which is specific to Lung cancer), so it seems the database is not only for GBM.
Question: I am using the following code to generate a tetramer from a monomer using a tetramer template in Modeller.

Template tetramer = 2jgs.pdb
Monomer = 2c4i.pdb

But my code only produces a monomer and not a tetramer.

env = environ()

#read pdb1 and pdb2 and creat and alignment between their sequences
aln = alignment(env)
mdl = model(env, file='2jgs.pdb')
aln.append_model(mdl=mdl, align_codes='2jgs', atom_files='2jgs')
mdl2 = model(env, file='2c4i.pdb')

aln.append_model(mdl=mdl2, align_codes='model1', atom_files='model1')
aln.write(file='2jgs_2jgs.ali')

z = (mdl.chains)
for i in z:
print i.name
print len(i.residues)

# Genereate and save structure=structure alignment
atmsel = selection(mdl.chains[i.name]).only_atom_types('CA')
r = atmsel.superpose(mdl2, aln,)
rms = r.rms
drms = r.drms

mdl2.write(file='model1_fit22.pdb')

Answer: I found the solution to be very simple. The mdl2.write() method actually overwrites the previous structures, so all I had to do was change the variable name of the output and join the files together
Question: Hello,

I have what might be a simple question but I would really like to get this right and despite all my readings it is still not clear to me what the best option is.

I am using DESeq2 to analyse an RNAseq experiment. The design of the experiment is as follow (each sample is in triplicate):

- cell line A - untrt
- cell line A - trt
- cell line B - untrt
- cell line B - trt

The goal is to peform the following contrats:

- cell line A - untrt **versus** cell line B - untrt
- cell line A - trt **versus** cell line B - trt
- cell line A - untrt **versus** cell line A - trt
- cell line B - untrt **versus** cell line B - trt

My question is, in order to do things "correctly", should I create a separate DESeq object for each contrast using only the samples I will be using for each contrast, or should I have just a single object from which I do the contrasts that I need?

For your information, the difference between these cell lines is that "cell line A" is overexpressing gene X, while "cell line B" is overexpressing genes X and Y.

Any help would be very much appreciated...

Answer: To make the comparisons you want, your sample data should include a column like cellline.treatment, with values like A.treat, A.untreat.

You make one dds object with cellline.treatment as the design with all the samples, and call results 4 times, each time specifying which cellline.treatment to compare to which.