cannin/biostars_qa huggingface dataset:
Question: Hi everyone,

For an allele mining project I need to generate in silico NGS data. I have chosen to use Grinder
to do so. I would like to add randomly generated SNPs to my reference genome, is any software
available to do that ?

Thanks for help!

Answer: Hi, you can also try Simulome, it is quite customisable.

https://github.com/price0416/Simulome

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5870732/
Question: I have a multiple samples and I run mapping with Bowtie2 for them and now I have .bam files. I have tried to run Multiqc to see how is mapping rate, duplication rate and other features but It did not work with bowtie output I think.
I would like to know is there any graphical program like Multiqc for observing aligning features?

Answer: MultiQC doesn't extract information from the bam files, rather, it extracts from log files. Bowtie2 logs are output to stderr, you have to redirect them to a file and pass this file to MultiQC.
Question: Is there any database or tool to get the list of pathways which a gene(*Homo sapiens*) is involved?
Like, Gene Symbol and list of pathways it is involved.

Answer: http://www.genecards.org/
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: https://gist.github.com/mdshw5/c7cf7a232b27de0d4b31
Question: I would like to know, can someone post a working example of converting a BAM file to CRAM and BAM again?

I tried samtools, cramtools but each time, the BAM tags get mangled. Can someone show me how to get:

samtools view file1.bam | md5sum
ed071dff8d60c74eefd1d14eb1f4bb21

then convert to CRAM, back to bam file2.bam such that:

samtools view file2.bam | md5sum
ed071dff8d60c74eefd1d14eb1f4bb21

Sorry if this is trivial but I have spent a few hours on this.

Answer: I think the source of the problem is that the BAM specification does not require or prescribe an order to the tags, hence these can be reordered and still represent the same information. I would contact the SAM developers and see what could be done about this:

https://github.com/samtools/
Question: Merry Christmas and Happy New Year,

Hi,I am new to this research.I am trying to understand how to fetch intron location.I have get gene location and exon location from NCBI, but I am still confused about exon location. My question could be shown as follows:

Structure1: **5'UTR -> exon -> intron ->exon->intron->exon-> 3'UTR**

Structure2: **5'UTR-> intron->exon->intron->exon->intron->3'UTR**

Q1: when I only consider intron location, whether I do not consider 5' or 3' UTR as intron? What I mean is 5'UTR or 3'UTR is not the same as intron.

Q2:As we learned, Structure1 is usually structure for gene.However, is that possible in Structure2? 5'UTR is adjacent with intron, but not exon.From NCBI .gff download file, I only could get gene location and exon location. If Structure2 is possible, how do I define intron location without 5'UTR location?

Thanks for your help.

Answer: Here we go, this is the way I see the definition a few of the terms used in your questions:

Exon: genomic sequence that is transcribed to ~~pre-mRNA~~ mature mRNA

Intron: genomic sequence that is not transcribed to ~~pre-mRNA~~ mature mRNA

UTR: mRNA sequence that is not translated to protein (hence, the name - un-translated region)

UTRs can be as short as a few bases to as long as multiple exons. (IIRC A4GALT1 has ~3 exons worth of 5' UTR - see [here][1])

So, we can conclude that a UTR can span an exon and hence can be followed by an intron. Both Structure 1 and Structure 2 are valid. Let's now tackle the problem once you're clear with this understanding. Exon co-ordinates should now give you UTR co-ordinates as well, but the challenge is to find out which part of the exons are UTR (AKA the translation start site position).

[1]: http://genome.ucsc.edu/cgi-bin/hgc?hgsid=402070773_F2iljTeLeSAU16aiC1Xcp2NRZaug&g=htcCdnaAliInWindow&i=NM_017436&c=chr22&l=43088126&r=43116876&o=43088126&aliTable=refSeqAli&table=refGene
Question: Hi,

The following piece of code uses Biopython to split all models from a valid PDB file into individual files.

It takes a PDB-formatted filename and a PDB code as input (arguments -f and -p respectively) and writes each existing model to a different output PDB file (with filename format "PDBcode"-"model".pdb).

My concern is that BioPython always numbers models starting from a dummy 0, so it is not directly possible to use the real model numbers in the name of the splitted files. (We can always add 1 to these numbers, but what if the PDB does not start from model 1, or model numbers are not consecutive?).

Is there any BioPython-compatible way of getting the actual model numbers from a PDB file and use these as part of the filenames for the spllited individual models?

Thanks!

#!/usr/bin/python
import sys, getopt
from Bio.PDB import PDBParser, PDBIO

def main(argv):
# Parse input arguments
pdbfile = ''
pdbid = ''
try:
opts, args = getopt.getopt(argv,"hf:p:",["pdbfile=","pdbid="])
except getopt.GetoptError:
print 'split_pdb_models.py -f -p '
sys.exit(2)

# Process input arguments
for opt, arg in opts:
if opt == '-h':
print 'split_pdb_models.py -f -p '
sys.exit()
elif opt in ("-f", "--pdbfile"):
pdbfile = arg
elif opt in ("-p", "--pdbid"):
pdbid = arg

# Split models from PDB into individual files
pdb = PDBParser(QUIET=True)
structure = pdb.get_structure(pdbid,pdbfile)
for model in structure:
io.set_structure(model)
io.save(structure.get_id() + "-" + str(model.get_id()) + ".pdb")
print 'Input file is:', pdbfile

if __name__ == "__main__":
io = PDBIO()
main(sys.argv[1:])

Answer: Yes, you can use the property "model.serial_num" to get the model number that was actually in the file (with no off-by-one shift), though I'm not sure this is documented anywhere.

See https://github.com/biopython/biopython/blob/master/Bio/PDB/PDBParser.py#L326-L334 and https://github.com/biopython/biopython/blob/20cb4d67adec136dcbe9bf01b29d56d39cdaa513/Bio/PDB/StructureBuilder.py#L69-L78 for where this happens in the code. It can also be seen in https://github.com/biopython/biopython/blob/1a5a47029ec69f839c0507fad0991c1f1a6ccbf3/Bio/PDB/Model.py that there is no "getter" function for this, but you can still access it directly with "model.serial_num".
Question: Of course, as a bioinformatician, I am aware of many large-scale open-source bioinformatics datasets, such as

- The ENCODE consortium (www.encodeproject.org, RNA-Seq, ChIP-Seq and so on),
- The Roadmap Epigenomics consortium (www.roadmapepigenomics.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The IHEC consortium (www.ihec-epigenomes.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The TCGA/ICGC consortia (www.cancergenome.nih.gov, www.icgc.org, large-scale cancer data, DNA-Seq, RNA-Seq, etc.) and
- The LINCS consortium (www.lincscloud.org/l1000 , gene expression for more than a million of different perurbation experiments).

I am wondering, however, what other wonderful datasets, the are both large and open-source, are currently available. That might include things like RNA-Seq, Chip-Seq, Bisulfite-Seq, whole genome sequencing, WGAS, and many other assays (not necessarily NGS-related, though that is what I am mostly looking for).

Also things like the (neural) connectome of certain species (in any event large data) could be of interest.

There are quite some GEO datasets that at least partially fulfill these requirements, but most are simply having to few data samples in order to be interesting to me.

Your comments are greatly appreciated!

Answer: - [European Nucleotide Archive][1]
- [COSMIC][2]
- [Phytozome][3]
- [JGI][4]

[1]: http://www.ebi.ac.uk/ena/browse
[2]: http://cancer.sanger.ac.uk/cosmic
[3]: http://phytozome.jgi.doe.gov/pz/portal.html
[4]: http://genome.jgi.doe.gov/
Question: I have a .csv file as follows:

,TEST1,TEST2
ENSG00000197421,2,0
ENSG00000213753,0,2
ENSG00000168746,0,2
ENSG00000261824,3,0
ENSG00000128310,1,2
ENSG00000235091,9,4

In R, I import the file like this:

> d <- read.csv("my_file.csv", header=TRUE, row.names=1)
> d
TEST1 TEST2
ENSG00000197421 2 0
ENSG00000213753 0 2
ENSG00000168746 0 2
ENSG00000261824 3 0
ENSG00000128310 1 2
ENSG00000235091 9 4

Checking that I do not have any duplicates:

> rownames(d)
[1] "ENSG00000197421" "ENSG00000213753" "ENSG00000168746" "ENSG00000261824"
[5] "ENSG00000128310" "ENSG00000235091"
> colnames(d)
[1] "TEST1" "TEST2"
> any(duplicated(rownames(d)))
[1] FALSE
> any(duplicated(colnames(d)))
[1] FALSE

Load libraries:

> suppressMessages(library("AnnotationDbi"))
> suppressMessages(library("org.Hs.eg.db"))

Then try to convert my Ensembl row names to Symbol in place:

> rownames(d) <- mapIds(org.Hs.eg.db,keys=rownames(d),column="SYMBOL",keytype="ENSEMBL",multiVals="first")
Error in `row.names<-.data.frame`(`*tmp*`, value = value) :
missing values in 'row.names' are not allowed

NOTE: Removing the first ',' on 'my_file.csv' did not help neither.

I managed to create a new field with the converted IDs but cannot replace it to the row names:

> d$SYMBOL <- mapIds(org.Hs.eg.db,keys=rownames(d),column="SYMBOL",keytype="ENSEMBL",multiVals="first")
> d
TEST1 TEST2 SYMBOL
ENSG00000197421 2 0 GGT3P
ENSG00000213753 0 2 CENPBD1P1
ENSG00000168746 0 2 LINC01620
ENSG00000261824 3 0 LINC00662
ENSG00000128310 1 2 GALR3
ENSG00000235091 9 4
> d_subset <- subset(d, !is.na(d$SYMBOL))
> d_subset
TEST1 TEST2 SYMBOL
ENSG00000197421 2 0 GGT3P
ENSG00000213753 0 2 CENPBD1P1
ENSG00000168746 0 2 LINC01620
ENSG00000261824 3 0 LINC00662
ENSG00000128310 1 2 GALR3
> rownames(d) <- d$SYMBOL
Error in `row.names<-.data.frame`(`*tmp*`, value = value) :
missing values in 'row.names' are not allowed

I don't get it.

Answer: In your last line, you're still trying to assign `rownames(d) <- d$SYMBOL`, and not your new `d_subset`.

`rownames(d_subset) <- d_subset$SYMBOL` should work.

Alternatively, if you wish to keep all the entries, you could retain the Ensembl names when it is not mapped:

> rownames(d) <- ifelse(is.na(d$SYMBOL), rownames(d), d$SYMBOL)

> d
TEST1 TEST2 SYMBOL
GGT3P 2 0 GGT3P
CENPBD1P1 0 2 CENPBD1P1
LINC01620 0 2 LINC01620
LINC00662 3 0 LINC00662
GALR3 1 2 GALR3
ENSG00000235091 9 4
Question: Dear all, I would like to know if there is any statistical test / program that allows me check whether the third codon position in my alignment is evolving too fast (too divergent) that it should not be included in a phylogenetic analysis. I googled around and didn't get any hits. I appreciate if you can share what you know on this topic.

Answer: There are a couple of tests for this, including Xia (2003, http://dx.doi.org/10.1016/S1055-7903(02)00326-3), which I gather is implemented in DAMBE.

Probably the most common way to asses this though, is to plot the number of 3rd codon substitutions in a sequence against some corrected distance . If the curve "plateaus" you know you are hitting saturation.
Question: Hello everyone,

Here's my problem: I have multiple VCF outputs that I have obtained using MuTect2 for 18 patients in total from Whole Exome Sequencing (WES/WXS) (so 18 VCF files). I have a BED file which has the list of all genes. FYI, MuTect2 reports single variants and INDELS in the VCF.

So, if I use *bedtools intersect*, it report each feature in A that overlaps B, even if its only **once**. What I want is to do an intersection between the list of genes in BED and all the VCF outputs, for example, at least 5 VCF that reports the same gene in the region. Bedtools describe here:

![intersectBed][1]

if I use:

***intersectBed -wa -a [BED] -b [all VCFs]***

It will report the list of all genes where the variant has been seen at least **once**, for example:

![intersect][2]

In the previous figure, you can see that bedtools (BT) will report gene B but I do not want that. What I want is that bedtools report gene A only if 5 VCF (or more) that have variants/INDELS overlap the same gene.

Is there any way to do it ? If there are other tools that you can suggest also, please let me know. Thank you in advance !

Alaa

[1]: http://bedtools.readthedocs.io/en/latest/_images/intersect-glyph.png
[2]: https://preview.ibb.co/drmtbQ/intersect.jpg

Answer: [BEDOPS `bedmap` and `vcf2bed`][1] will do what you want:

$ bedmap --count --echo --echo-map-id-uniq --delim '\t' genes.bed <(vcf2bed < snps.vcf) | awk '($1>=5)' | cut -f2- > answer.bed

You can further constrain the VCF conversion step with `--snvs`, `--insertions`, and `--deletions`. See the documentation in `vcf2bed --help` or [online][2] for more detail.

If you have a lot of VCF files and you want to collate them into one file, so that you can map them in one pass, you can union them into a single BED with `bedops --everything`:

$ bedops --everything <(vcf2bed < A.vcf) <(vcf2bed < B.vcf) ... <(vcf2bed < N.vcf) > unionSNPs.bed

Then you would map against the unioned SNPs file:

$ bedmap --count --echo --echo-map-id-uniq --delim '\t' genes.bed unionSNPs.bed | awk '($1>=5)' | cut -f2- > answer.bed

Either way, putting `--count` first in the list of options to `bedmap` prints the number of overlapping SNPs in the first column of the output. The `awk` step filters on this value, printing any mapped elements that have 5 or more overlapping SNPs. Then the `cut` step removes the column with the number of overlaps, leaving a proper BED file as output.

The file `answer.bed` will have the gene and the IDs of any overlapping SNPs. If you want the SNP information in entirety, or their positions, you can use other `--echo-map-*` options. Or if you don't care about the SNPs, you can leave out `--echo-map-id-uniq` and just print the gene with `--echo`. See the documentation in `bedmap --help` or [online][3] for more detail.

[1]: https://github.com/bedops/bedops
[2]: http://bedops.readthedocs.io/en/latest/
[3]: http://bedops.readthedocs.io/en/latest/
Question: Hi all,
I've just written a tool adding one or more **extra column** in a [VCF file](http://www.1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcf3.3). The header now looks like this:

(...)
#CHROM POS ID REF ALT QUAL FILTER INFO MY_COL1 MY_COL2 FORMAT NA00001 NA00002 NA00003
(...)

Is there something in the **VCF** spec saying that another column can't be added ?
because when I used **VCFTOOLS**, it says:

vcftools --vcf file.vcf
(...)
Scanning file.vcf ...
Ninth Header entry should be FORMAT: MY_COL1
Currently scanning CHROM: 19
Currently scanning CHROM: 20
Currently scanning CHROM: X

Answer: Which kind of information do you want to add? I don't think that the VCF specifications allow to add new columns, but since this format is still in an early phase of development, you could contact the authors and propose them a new functionality.

However, VCF files should be used only to describe the SNPs and their genotypes, and any other kind of information should go somewhere else... for example, if you have statistics associated with a snp, you should consider a flat file or a database.
Question: I notice GSEA automatically geneartes a heat-map like diagram for the pathways but does not include a 'scale' bar.
I tried to get a better understand on how the heat map is plotted, but the explanation is not very clear on the tutorial site.
Does anyone know how is the heat map generated? Are the colors based on fold change, or something else?

Thank you

Answer: Hi! My understanding is that expression levels are visualized by colors, with high expression shown in red and low expression shown in blue:
"In a heat map, expression values are represented as colors, where the range of colors (red, pink, light blue, dark blue) shows the range of expression values (high, moderate, low, lowest)."
Question: Hello.

I read the paper "Whole-exome sequencing identifies a recurrent NAB2-STAT6 fusion in solitary fibrous tumors" briefly.

After reading it, I tried to find the gene fusion in my exome sequence data with cancer patients from hospital.

But I think that it is a better strategy to identify the gene fusion case through the use of the public data that is specified in the above paper before trying to test on my own set.

So I tried to find an explanation about how to access public data in the paper to identify the gene fusion.

Finally I found the bold sentence (below).

**Accession code. Sequence data used for this analysis are available at the database of Genotypes and Phenotpyes(dbGap) under accesssion. phs000568.v1.p1.**

Just reading it, I think I can download data if i have permission. So I go to dbGap stite and search phs000568.v1.p1 and many results occur to me. But I don't know how to access above databases and also don't know how to download.

So Is there anyone who have experienced this before?

If so, I need your help!

Thank you.

Answer: You can find the relevant [dbgap][1] study record by [searching dbgap for accession "phs000568.v1.p1"][2]. This will get you to the appropriate dbgap study page:
http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000568.v1.p1

From there we can learn that the study is available for "General Research Use" and does not require IRB. This makes the process much easier than it would be otherwise. But, you still need to have a PI apply for permission to access the data and they will have to agree to follow/supervise the proper protocols in terms of how the data is used and stored.

If your PI is USA-based and has applied for funding then they will most likely have an ERA commons account and be able to login and initiate request for access. You will then probably have to wait some time for your institution signing official and then dbgap to approve your request. I have no idea how this works for non-USA based researchers though.

[1]: http://www.ncbi.nlm.nih.gov/gap
[2]: http://www.ncbi.nlm.nih.gov/gap/?term=phs000568.v1.p1
Question: Hi,

I have a fasta file **input.fa** which has duplicate fasta headers.

>UID584_Org_name_strain
JGTTQEWEILYNMCCDASSHHGTQEFGRKKQLAADBNGTAVQEGFN
>JGH1236_Org_name
HFHGTSNNIGLTTREKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLVDGNVLTDRLSIGGKTALTGVDP
>JGH1236_Org_name
FHGTSNNIGLTTREKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLVDGNVLTDRLSIGGKTALTGV
>KIL563.2_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSNNIGLTTREKKQLAADBN
>KIL563.2_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLV
>KIL563.2_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGFNMENDLALFVTY
>GTK584_Org_name_str
KKQLAADBNGTAVQEGFNMENDLALFVTYGHFHGTSNNIGLTTREKKQLAAWLMAFDSCVIIPPTRE
>GTK584_Org_name_str
DBNGTAVQEGFNMENDLALFVTYGHFHGTSNNIGLTTREKKQLAAWLMAFDSCVIIPPTREADBNGTAVQEGFNME
>EAD5624_Org_nam
LTTREKKQLAADBNGTAVQEGFNMENDLALFVTYBNGTAVQEGFNMENDLALFVTYGAAWLMAFDSCVIIPPT
>EAD5624_Org_nam
LTTREKKQLAADBNGTAVQEGFNMENDLALFVTYBNGTAVQVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGF

I want to rename all the headers with numbers before the first underscore '_' like this:

>UID584.1_Org_name_strain
JGTTQEWEILYNMCCDASSHHGTQEFGRKKQLAADBNGTAVQEGFN
>JGH1236.1_Org_name
HFHGTSNNIGLTTREKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLVDGNVLTDRLSIGGKTALTGVDP
>JGH1236.2_Org_name
FHGTSNNIGLTTREKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLVDGNVLTDRLSIGGKTALTGV
>KIL563.2.1_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSNNIGLTTREKKQLAADBN
>KIL563.2.2_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGFNMENDLALFVTYGAHLV
>KIL563.2.3_Org_name
TTFWLMAFDSCVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGFNMENDLALFVTY
>GTK584.1_Org_name_str
KKQLAADBNGTAVQEGFNMENDLALFVTYGHFHGTSNNIGLTTREKKQLAAWLMAFDSCVIIPPTRE
>GTK584.2_Org_name_str
DBNGTAVQEGFNMENDLALFVTYGHFHGTSNNIGLTTREKKQLAAWLMAFDSCVIIPPTREADBNGTAVQEGFNME
>EAD5624.1_Org_nam
LTTREKKQLAADBNGTAVQEGFNMENDLALFVTYBNGTAVQEGFNMENDLALFVTYGAAWLMAFDSCVIIPPT
>EAD5624.2_Org_nam
LTTREKKQLAADBNGTAVQEGFNMENDLALFVTYBNGTAVQVIIPPTREWWQQLGTGTSEKKQLAADBNGTAVQEGF
It would be convenient if anybody could suggest me `sed/awk/grep` command for this.
Any help would be appreciated. Thanks!

Answer: This will fail if the sequences are not ordered as in your example:

awk 'BEGIN{OFS=FS="_"}{if(/^>/){CUR=$1;{if(CUR==PRE){NUM++}else{NUM=1}};$1="";print CUR"."NUM $0;PRE=CUR}else{print $0}}' in > out
Question: I have two wig files and wish to perform correlation plot before merging the two files. I checked with http://cistrome.org/ for performing the needful. Was interested to know more ways to perform it.

Answer: You can use the `wigToBigWig` and `bigWigCorrelate` tools from UCSC (find the executable for your platform [here][1]).

First convert your WIG files to bigWig:

```
wigToBigWig v 4 - Convert ascii format wig file (in fixedStep, variableStep
or bedGraph format) to binary big wig format.
usage:
wigToBigWig in.wig chrom.sizes out.bw
Where in.wig is in one of the ascii wiggle formats, but not including track lines
and chrom.sizes is two column:
and out.bw is the output indexed big wig file.
Use the script: fetchChromSizes to obtain the actual chrom.sizes information
from UCSC, please do not make up a chrom sizes from your own information.
options:
-blockSize=N - Number of items to bundle in r-tree. Default 256
-itemsPerSlot=N - Number of data points bundled at lowest level. Default 1024
-clip - If set just issue warning messages rather than dying if wig
file contains items off end of chromosome.
-unc - If set, do not use compression.
-fixedSummaries - If set, use a predefined sequence of summary levels.
-keepAllChromosomes - If set, store all chromosomes in b-tree.
```

Then correlate them:

```
​bigWigCorrelate - Correlate bigWig files, optionally only on target regions.
usage:
bigWigCorrelate a.bigWig b.bigWig
options:
-restrict=restrict.bigBed - restrict correlation to parts covered by this file
-threshold=N.N - clip values to this threshold
```

[1]: http://hgdownload.cse.ucsc.edu/admin/exe/