cannin/biostars_qa huggingface dataset:
Question: I'm trying to write a CWL workflow that will do a simple trim and map procedure on some bacterial short read data that I've got. For the trimming I'm planning on using trimmomatic and I have modified the Duke-GCB trimmomatic tool slightly for this purpose (see https://github.com/pvanheus/GGR-cwl/blob/master/trimmomatic/trimmomatic.cwl). Here is the CWL for the workflow (`trim.cwl`):

cwlVersion: v1.0
class: Workflow
inputs:
reads1: File
reads2: File
slidingw: string
minl: string
outputs:
trimmed_reads:
type: File
outputSource: trim_reads/output_read1_trimmed_file

steps:
trim_reads:
run: { $import: "../../tools/trimmomatic/trimmomatic.cwl" }
inputs:
input_read1_fastq_file: reads1
input_read2_fastq_file: reads2
slidingwindow: slidingw
minlen: minl
out: [output_read1_trimmed_file]

but trying to run this with `cwltool trim.cwl` yields:

cwltool 1.0.20160714182449
Tool definition failed validation:
Validation error in object file:///net/ceph-mon1.sanbi.ac.za/sanbi/scratch/pvh/cwl/workflows/trim/trim.cwl
Could not validate as `CommandLineTool` because
could not validate field `outputs` because
At position 0
could not validate field `outputSource` because it is not recognized and strict is True, valid fields are: label, secondaryFiles, format, streamable, doc, id, outputBinding, type

missing required field `baseCommand`
could not validate field `steps` because it is not recognized and strict is True, valid fields are: id, inputs, outputs, requirements, hints, label, doc, cwlVersion, class, baseCommand, arguments, stdin, stderr, stdout, successCodes, temporaryFailCodes, permanentFailCodes
Could not validate as `ExpressionTool` because
could not validate field `outputs` because
At position 0
could not validate field `outputSource` because it is not recognized and strict is True, valid fields are: label, secondaryFiles, format, streamable, doc, id, outputBinding, type

missing required field `expression`
could not validate field `steps` because it is not recognized and strict is True, valid fields are: id, inputs, outputs, requirements, hints, label, doc, cwlVersion, class, expression
Could not validate as `Workflow` because
could not validate field `steps` because
the value `[{'id': 'file:///net/ceph-mon1.sanbi.ac.za/sanbi/scratch/pvh/cwl/workflows/trim/trim.cwl#trim_reads',
'inputs': [{'id': u'file:///net/ceph-mon1.sanbi.ac.za/sa[...]`
is not a valid type in the union, expected one of:
- array of , but
At position 0
missing required field `in`
could not validate field `inputs` because it is not recognized and strict is True, valid fields are: id, in, out, requirements, hints, label, doc, run, scatter, scatterMethod

I'm sure there are all sorts of errors in my CWL, but the parsing as a CommandLineTool is surely the wrong thing?

Answer: Hello Peter vH!

For workflows we have "in" and "out". CommandLineTools use 'inputs' and 'outputs'. You've mixed the two in your workflow and the last line of the overly detailed error message says as much. Our apologies for the error messages being so hard to read.
Question: Hi all,

Sorry I know that this question has been asked several times, but unfortunately I haven't been able to find the right answer, or didn't understand.

I'm trying to get TMM normalized counts thanks to edgeR.

I understand that I have to compute normalization factors :

dgList <- calcNormFactors(dgList, method="TMM")

which gives me a normalization factor for all samples :

head(dgList$samples)

group lib.size norm.factors
S1 1 21087314 0.9654794
S2 1 16542810 1.1589117
S3 1 18875473 0.8763291
S4 1 15865414 1.0864038
S5 1 19179795 1.0488230
S6 1 15063992 1.0707007

But at this step I don't know what to do to get a matrix of normalized TMM counts.

I know that I can get CPM normalized counts thanks to :

cpm(dgList)

But CPM and TMM are not the same, right ?

Thanks in advance for any of your input on this topic.

Answer: If you run the cpm function on a DGEList object which contains TMM normalisation factors then you will get TMM normalised counts. Here is a snippet of the source code for the cpm function:

cpm.DGEList <- function(y, normalized.lib.sizes=TRUE, log=FALSE, prior.count=0.25, ...)
#	Counts per million for a DGEList
#	Davis McCarthy and Gordon Smyth.
#	Created 20 June 2011. Last modified 10 July 2017
{
lib.size <- y$samples$lib.size
if(normalized.lib.sizes) lib.size <- lib.size*y$samples$norm.factors
cpm.default(y$counts,lib.size=lib.size,log=log,prior.count=prior.count)
}

The function checks to see if a DGEList object was provided with a *lib.size* and *norm.factors* column (created when you run calcNormFactors), if so then it uses those in the normalisation of the raw counts. You were right in your original post, just run the following and you will have TMM normalised counts:

dge <- calcNormFactors(dge, method = "TMM")
tmm <- cpm(dge)
Question: I've simulated RNA abundance with wgsim. The simulation itself was error free. There is a single factor in my experiment that looks like:

```r
A1 A2 A3 B1 B2 B3
R1_101 113 113 113 13 11 9
R1_102 247 246 246 12 12 14
R1_103 20835 20915 20788 9973 9955 9973
```

A1, A2, A3 are the simulated replicates for the first level. B1, B2 and B3 are the simulated replicated for the second level.

As expected, the reads counts for each level are very close because it was an error-free simulation. The purpose of the experiment is to compare it with cuffdiffs (another differential package) in detecting log-fold changes.

Unfortunately, I run into an error in DESq2:

```r
Error in estimateDispersionsFit(object, fitType = fitType, quiet = quiet) :
all gene-wise dispersion estimates are within 2 orders of magnitude
```

It looks like the package is unable to estimate a dispersion factor (most likely it's too small). However, I had no problem with cuffdiffs. Is there anything that I can do to make it work?

Answer: The DESeq2 method has empirical Bayes parts to it, which involve sharing information across genes to improve estimates (see [the paper][1]). In this case, during dispersion estimation, we look across the genes at the distribution of gene-wise dispersion estimates to improve the final estimates (posterior modes). If you simulate data which has no (over)dispersion, then these methods don't make sense.

Did you clip the warning message which was printed in the console for some reason? Because the rest of it tells you what to do:

```r
# all gene-wise dispersion estimates are within 2 orders of magnitude
# from the minimum value, and so the standard curve fitting techniques will not work.
# One can instead use the gene-wise estimates as final estimates:
dds <- estimateDispersionsGeneEst(dds)
dispersions(dds) <- mcols(dds)$dispGeneEst
```

[1]: http://www.genomebiology.com/2014/15/12/550
Question: I would like to identify a "canonical" transcript for every protein-coding gene in Ensembl. For project-related reasons, I'm using the `EnsDb.Hsapiens.v75` package in R. I realize, of course, that "canonical" is a working definition at best, and inappropriate in some cases - but for ease of graphing some data I just want one transcript per gene for now. From manually inspecting genes in Ensembl, it looks like the lowest-numbered transcript ID for each corresponds to what I'm looking for. Some code to pull out a few examples:

library(EnsDb.Hsapiens.v75)
library(tidyverse)

genes <- keys(EnsDb.Hsapiens.v75, keytype='GENEID')
ensembl <- AnnotationDbi::select(EnsDb.Hsapiens.v75, keys=genes, keytype='GENEID',
columns=c('GENEID', 'SYMBOL', 'GENEBIOTYPE'))
ensembl_cds <- filter(ensembl, GENEBIOTYPE=='protein_coding')
ensembl_cds_tx <- AnnotationDbi::select(EnsDb.Hsapiens.v75, keys=genes, keytype='GENEID',
columns=c('SYMBOL', 'TXID'))
head(ensembl_cds_tx)

gois <- c('RSPO1', 'PRSS1', 'CDH1')
gois_tx <- filter(ensembl_cds_tx, SYMBOL %in% gois) %>% arrange(SYMBOL, TXID) %>% print()
gois_tx_lowest <- gois_tx[!duplicated(gois_tx$SYMBOL),] %>% print()
Each of the lowest transcript IDs pulled out above (ENST00000261769, ENST00000311737, ENST00000356545) corresponds to an Ensembl transcript for the respective genes (CDH1, PRSS1, RSPO1) that matches with Refseq and the Consensus CDS database. (Although, for RSPO1, there are three other transcripts that also have Refseq matches, which speaks to the arbitrariness of picking a single canonical transcript.)

My question is, is this the general practice across the Ensembl transcript database, that the lowest numbered transcript for a gene corresponds to a canonical or semi-canonical transcript, or have I just gotten lucky so far?

Answer: No. The numbers are arbitrary. The canonical transcript is the one which is labelled canonical, which you can get as a filter or an attribute.

The stable IDs are assigned in order, so the first transcript every identified was ENST00000000001, the second ENST00000000002 etc. This means that for a gene, the one with the lowest number was the first one to be identified. In all probability, the first one identified is the one that is the most highly expressed, highly conserved and well-studied, which makes it coincidentally also the canonical. But it's not always the case.
Question: Hi, I am trying to do PCA analysis for my samples for initial quality control. I have 2 different sets of samples - one were sequenced 50bp and other was sequenced 75bp (both of them have disease and control cases). To do the PCA on those samples, I ran DEseq2 on them (which necessarily requires non-normalized counts), followed by vst and plotPCA. But in my PCA plot, I see two clusters - one for the 50bp samples and the other for 75 bp samples. This is not necessarily expected, since there is nothing different between the samples except for sequencing depth. Someone said I should normalize my data. But I think that will be taken care of by DESeq2, and it anyway shouldn't be fed normalized counts. Here is the plot.

Any suggestions?
![PCA plot link][1]

Here is the image link in case it doesn't show - https://imgur.com/CFXHky4

Here is my code. Here P and NP mean Pain and Non-Pain, which is the effect I am studying. Some files are 50bp, others are 75bp, and I have **not** provided that info to DESeq -

counts_all_fullGTF <- featureCounts(nthreads=3, isGTFAnnotationFile=TRUE, annot.ext='/Volumes/bam/DRG/annotations/Homo_sapiens.GRCh38.95.gtf', files=c('/Volumes/.../47T7L.fastqAligned.sortedByCoord.out.bam','/Volumes/.../47T7R.fastqAligned.sortedByCoord.out.bam'))$counts

sampleTable_all <- data.frame(condition=factor(c('P','NP','P','P','P','P','NP','P','P','NP','P','P','P','P','P','P','NP','P','P','NP','P','P','P','P','P','NP','NP','NP','NP','P','P','P','P','P')))

coldata <- sampleTable_all
deseqdata_fullGTF <- DESeqDataSetFromMatrix(countData=counts_all_fullGTF, colData=coldata, design=~condition)
dds_fullGTF <- DESeq(deseqdata_fullGTF)
vsd_fullGTF <- vst(dds_fullGTF)
library(ggrepel)
plotPCA(vsd_fullGTF, ntop=5000, "condition") + geom_label_repel(aes(label=substr(name, start = 1, stop = 6), colour = "white", fontface = "bold"))

[1]: https://imgur.com/CFXHky4

Answer: Firstly, yes, the data should be normalised, and ideally transformed to a normal distribution.

When you ran `vst()`, did you set `blind = FALSE` and did you have `ReadLength` (or `batch`) as a covariate in your design formula? Setting `blind = FALSE` 'exposes' the variance stabilising transformation to the design formula - it's useful for producing expression data for downstream applications outside of DESeq2. Read here for further information: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#blind-dispersion-estimation

If you still need to remove that batch effect, then use `limma::removeBatchEffect()` on the variance-stabilised expression levels.

As to the origin of the observed difference, I think that multi-mapping will obviously be a bigger issue, given the shorter read length. However, multi-mapping may be dealt with by the aligner or pseudo-aligner that you use - which were used?

Kevin
Question: `**Package ‘![VennDiagram][1]’.**` How can I get the list of common objects in each combination?

For e.g. what are the common 44 Genes in interset (similar image got from online)?

I have used below scripts

Overlap <- calculate.overlap(x = list(BB, BB, BRR))

Overlap[[1]]

But in this method, I can not find the list of 1761, 466 genes

how these numbers are given?

Even when I tried `Overlap$a5` , `Overlap$a7`, etc this also not given some of the values correctly?

Could someone help me?

Answer: From what I can tell, the `a`s are distributed from top left to right, to bottom.
i.e.
```
a1 = 1761
a2=126
a3=466
a4=64
a5=44
a6=27
a7=366
```

This is based on:
![enter image description here][1]

> library("VennDiagram")
>
> gene_list = paste0("GENE", 1:1000)
>
> studies = list( S1=sample(gene_list, 700, replace = FALSE),
+ S2=sample(gene_list, 700, replace = FALSE),
+ S3=sample(gene_list, 700, replace = FALSE) )
>
> ol = calculate.overlap(x = studies)
> ol_size=sapply(ol, length)
>
>
> venn.diagram(
+ x = studies,
+ euler.d = TRUE,
+ filename = "Euler_3set_scaled.tiff",
+ cex = 2.5,
+ cat.cex = 2.5,
+ cat.pos = 0
+ );
[1] 1
> # Get the 63
> length( setdiff(studies$S1, union(studies$S2, studies$S3)) )
[1] 63
> # Confirm order
> ol_size
a5 a2 a4 a6 a1 a3 a7
351 143 143 134 63 72 72

And it matches the order from `a1` to `a7` going left-right, top-bottom.

If you want to be sure of what you're getting, you can use [set operators][2] like I did with `setdiff()` and `union()`.

[1]: https://i.imgur.com/85fPa1e.png
[2]: https://stat.ethz.ch/R-manual/R-devel/library/base/html/sets.html
Question: I'm in the process of creating concensus sequences, and my specific aim is to remove any ambiguous characters (for reasons that are unimportant to the question).

I couldn't find a simple tool that would do what I want satisfactorily, so I started writing my own (until someone reminded me about HMMs but more on that is a second).

I have the following MSA as a AlignIO object:

>>> print(msa)
SingleLetterAlphabet() alignment with 16 rows and 149 columns
MSTTPEQIAVEYPIPTYRFVVSLGDEQIPFNSVSGLDISHDVIE...QAA PAU_02775
MSTTPEQIAVEYPIPTYRFVVSIGDEQIPFNSVSGLDISHDVIE...QAA PLT_01696
MSTTPEQIAVEYPIPTYRFVVSIGDEQVPFNSVSGLDISHDVIE...QAA PAK_02606
MSTTPEQIAVEYPIPTYRFVVSIGDEKVPFNSVSGLDISHDVIE...QAA PLT_01736
MTTTT----VDYPIPAYRFVVSVGDEQIPFNNVSGLDITYDVIE...QAA PAK_01896
MATTT----VDYPIPAYRFVVSVGDEQIPFNSVSGLDITYDVIE...QAA PAU_02074
MSVTTEQIAVDYPIPTYRFVVSVGDEQIPFNNVSGLDITYDVIE...QAA PLT_02424
MTITPEQIAVDYPIPAYRFVVSVGDEKIPFNNVSGLDVHYDVIE...QAP PLT_01716
MAITPEQIAVEYPIPTYRFVVSVGDEQIPFNNVSGLDVHYDVIE...QAA PLT_01758
MSTSTSQIAVEYPIPVYRFIVSIGDDQIPFNSVSGLDINYDTIE...QAV PAK_03203
MSTSTSQIAVEYPIPVYRFIVSVGDEKIPFNSVSGLDISYDTIE...QAV PAU_03392
MSITQEQIAAEYPIPSYRFMVSIGDVQVPFNSVSGLDRKYEVIE...QVP PAK_02014
MSITQEQIAAEYPIPSYRFMVSIGDVQVPFNSVSGLDRKYEVIE...QVP PAU_02206
MSTTADQIAVQYPIPTYRFVVTIGDEQMCFQSVSGLDISYDTIE...EFH PAK_01787
MSTTADQIAVQYPIPTYRFVVTIGDEQMCFQSVSGLDISYDTIE...EFH PAU_01961
MSTTVDQIAVQYPIPTYRFVVTVGDEQMSFQSVSGLDISYDTIE...EFH PLT_02568
*

(Note the asterisk I've added for the moment).

If I build a concensus with `hmmemit -c` ("majority-rule consensus sequence"), I get the following sequence:

>hmmemit-consensus *
MSTTAEQIAVEYPIPTYRFVVSVGDEQIPFNSVSGLDISYDVIEYKDGVGNYYKMPGQRQ
AINITLRKGVFSGDTKLFDWINSIQLNQVEKKDISISLTNEAGTEILLTWSVANAFPTSL
TSPSFDATSNEVAVQEISLTADRVTIQAA

At column 23 (22 if zero based) (the asterisk), `hmmemit` places a Valine (`V`). My own tool however, places an Isoleucine (`L`).
Part of my process is to use `collections.Counter` to score the string, and in that column, `I` occurs 8 times, and `V` 7 times. Why is `hmmemit` choosing a lower frequency amino acid?

Basic process of my script (relevant parts only):

def enumerate_string(string):
"""Returns the most common characters of a string. Multiple characters are returned if there are equally frequent characters."""

from collections import Counter
counts = Counter(string)

keys = []
for key, value in counts.iteritems():
if value == max(counts.values()):
keys.append(key)

return keys

------

>>> msa = AlignIO.read('myalignment.fasta', 'fasta')
>>> msa_summary = AlignInfo.SummaryInfo(msa)

>>> Counter(msa_summary.get_column(22))
Counter({'I': 8, 'V': 7, 'L': 1})
>>> enumerate_string(msa_summary.get_column(22))
['I']

This clearly shows that Isoleucine is the most common, so why doesn't `hmmemit` choose it?

Answer: Because hmmemit is emitting the maximum likelihood sequence, according to the profile HMM's parameterization.

("majority-rule" in hmmemit means w.r.t. the profile HMM, not w.r.t. the MSA that the HMM was built from. hmmemit takes the HMM as input, not the MSA.)
Question: Dear all, there are many posts about remove duplicate sequences in a fasta file (https://www.biostars.org/p/3003/), but I want to remove only the duplicate sequences with the same ids.

I have many duplicate sequences in my fasta file, but with different ids and I want to keep them.

How to remove only same id sequence duplicates? I have protein sequences and my sequences are split in different lines.

Answer: http://bioinf.shenwei.me/seqkit/usage/#rmdup
Question: Hi everyone,

I have 2 datasets that I will be using for my university project. The first dataset is relatively small and the other is much larger. Would it look odd in my university project to have a different threshold for the interaction probability for each of the data sets. The first one is smaller and I set the interaction probability to 0.5. Whereas the second dataset has several lists of proteins that are nearly 1000 proteins in length, so I will have to up the interaction probability threshold to 0.7 or 0.9 to make my Cytoscape/STRING network look less 'messy'. Is this bad practice?

Answer: While it is usually best and easier to explain/justify the use of a single threshold, I think it is ok to use different ones as long as it is properly explained. The true bad practice here would be to use different threshold and not make it clear. One way to be completely transparent about this would be to build 4 networks: the big and small dataset with both the high and low threshold. Then it is easy to discuss the benefits of each threshold by showing both of them and by explaining, for instance, why the low threshold is better for visualization of the small dataset and the high threshold is better for the visualization of the big dataset.
Question: Hi all.

Now, I’m trying to get an assembly stat with reference to table S1 and figure S1 in this paper.
> [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7266049/][1]

I used the following command to get the result of BLAST for my fasta.

blastx -query /home/nkarim/avenae/trinity_even_out_dir/Trinity.300.longest.fasta \
-db /home/nkarim/blast/db/nr \
-outfmt 6 \
-evalue 1e-3 \
-out /home/nkarim/blast/output/avenae.out

But I don't know how do I get the items;"Top BLASTx-hit species", "Percent of gene with at least one BLASTx hit" and "Percent of gene with at least one GO annotation".
Can I get each items using the result of BLAST?

And could you tell me how to get a table like table S1 if you know it?

Thank you very much for your help in advance!

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7266049/

Answer: Since you're blasting with ncbi-NR you're almost there. By changing the blastx output options to outfmt 6 you're receiving tabular output in `avenae.out`, you can customise the output columns to also receive the species names.

For example, you could run:
```
blastx -query /home/nkarim/avenae/trinity_even_out_dir/Trinity.300.longest.fasta \
-db /home/nkarim/blast/db/nr \
-outfmt "6 qseqid sseqid sscinames scomnames staxid pident length mismatch gapopen qstart qend sstart send ppos evalue bitscore" \
-evalue 1e-3 \
-out /home/nkarim/blast/output/avenae.out
```

Now you have all these extra columns, `sscinames` and `scomnames` are useful for your end as these contain the scientific name and the common name of your hits.

The laziest way to get the 'top' hit of your queries is to keep only the first hit seen, that normally has the lowest e-value or the highest score :

awk '! a[$1]++' avenae.out > avenae_first_hit_only.out

This falls apart when you have several species with equally good hits, you should still check whether the first hits reported are actually the best hits (manually?). From this out file you also get the percentage of genes with at least one BLAST hit:

wc -l avenae_first_hit_only.out

That will report the number of genes with at least one hit. For GO terms, I'm not sure - you could run GO annotation for your own genes, there are a few online tools for that like PANNZER2 or eggnog-mapper where you can upload your proteins.

For Table S1 (Assembly stats), Trinity has scripts for that, see this older thread https://www.biostars.org/p/233160/
Question: I am trying to calculate the MFE along with the secondary structure for multifasta using RNAFold. The output generated is of the format.

>abc
GGCGGAGGUAGGGAGGCACGCGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAG
(((((.....(((((((.(..((.......)).).)))))))........((((((.((............)).)))))).((((....))))))))).. (-35.80)
>lmn
GGGAGGCACGCGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAGUACCACCCCA
(((((((.(..((.......)).).))))))).............((((..((((.........((((.(.((((....)))).).)))))))))))).. (-29.30)
>xyz
CGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAGUACCACCCCACCCCGGGACA
....(((........)))(((((............((((.(((((.((......((((.(.((((....)))).).)))))).))))).))))))))).. (-28.40)

Is there a way to get the output in tubular format with 1. Identifier, 2. sequence, 3. secondary structure and 4. MFE as columns?
I have written regular expression scripts to capture each of the four and paste it in a file but I don't think that's an efficient way of doing it. Is there any other convenient way of doing it?

Answer: In Perl:

#!/usr/bin/perl

use strict;
use warnings;

while (<>) {
chomp;
if (/>/) {
s/>//;
print "$_\t"; # just the seq id
}
elsif (/\((-\d+\.\d+)\)$/) {
my $mfe = $1;
s/ \($mfe\)//;
print "$_\t$mfe\n"; # fold + MFE
}
else {
print "$_\t"; # the seq
}
}

Validation:

$ perl fold2tab.pl < fold.fa
abc GGCGGAGGUAGGGAGGCACGCGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAG (((((.....(((((((.(..((.......)).).)))))))........((((((.((............)).)))))).((((....))))))))).. -35.80
lmn GGGAGGCACGCGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAGUACCACCCCA (((((((.(..((.......)).).))))))).............((((..((((.........((((.(.((((....)))).).)))))))))))).. -29.30
xyz CGAUGGUAUUUCAGAGCCUCCCGAAUACAACUCCAGGGUAGGGUGUUGAAAGCGUUGGAGAUGUCUAAAGACACCGCCAGUACCACCCCACCCCGGGACA ....(((........)))(((((............((((.(((((.((......((((.(.((((....)))).).)))))).))))).))))))))).. -28.40
Question: I'm a new PhD student and bioinformatics is all very new to me so apologies if this seems trivial. I have some paired-end RNA-seq data which I've managed to convert to sorted and indexed bam files however I'm confused to where I go from here, am I right in thinking I can now visualise this data and if so how?

Answer: You can visualise the indexed BAMs with [IGV][1]. It should be pretty straight forward if you **choose the correct genome build**.

If you are working human or mouse genome, or any of the genomes listed on [ucsc genome browser][2], you can convert your bam files to bigwig files using [deepTools2][3] and create [tracks][4] to be viewed on genome browser. These tracks can be saved and visualised on any web browser and can be shared with collaborators.

[1]: http://software.broadinstitute.org/software/igv/home
[2]: https://genome-euro.ucsc.edu/cgi-bin/hgGateway?hgsid=602554305_kMaKYkyWlt61jjMmoBrR3w3UtOHA&redirect=manual&source=genome.ucsc.edu
[3]: http://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html
[4]: https://genome.ucsc.edu/goldenpath/help/customTrack.html