# Question
Title: How can I calculate coverage at single bases using a bam file?
I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:

```
chr base_position total_coverage fwd_coverage rev_coverage

```

I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information.

# Answer
*Note: not yet tested, so there may be some additional fiddling with command line options needed*

The per-base depth can be obtained from *samtools depth* (-a includes zero-coverage positions):

```
samtools depth -a in1.bam > depth_in1_both.tsv

```

To split this by forward and reverse, you can use an initial pipe through *samtools view* to exclude or include reverse-complement mappings:

```
samtools view -b -F 0x10 in1.bam | samtools depth -a > depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a > depth_in1_rev.tsv

```

To filter for particular regions, this can be pre-piped through *bedtools intersect*:

```
bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a > depth_in1_regions_fwd.tsv

```

Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...

If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using `samtools depth -aa`), then the lines can be *paste*d, with anonymous pipes, with unnecessary fields filtered out using *cut*:

```
paste <(samtools depth -aa in1.bam) \
      <(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      <(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

```

Or equivalently, if using intermediate files that have been generated already:

```
paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

```

It would also be possible to use *join* to combine the files (in cases where file output is *not* consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.
> 4 votes

# Answer
```
git clone https://github.com/lh3/htsbox
cd htsbox && make
./htsbox pileup -cCf ref.fa aln.bam | less -S

```

This output a VCF containing positions covered by at least one read. Tag `ADF` and `ADR` give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.
> 4 votes

# Answer
Brent Pedersen (@brentp) claims mosdepth is 2x as fast as `samtools depth`. I'm not sure whether this provides strand information, but it might be the fastest tool available.

Another tool of possible interest (also written by Brent) is indexcov. This doesn't give single bp precision, but is able to provide even quicker estimates of coverage at 16Kb resolution.

Installing these tools from source is non-trivial, and I struggled with installation in the past. However, they are both now available in bioconda, which should make installation simple.
> 0 votes

---
Tags: sam, bedtools, coverage
---


# Question
Title: Can we use GTEx data as control data for TCGA data?
I am using Recount2 TCGA data and was wondering is it right to use GTEx data as control data for this. I would really appreciate your views on this?

# Answer
Definately not directly. Probably not even if used indirectly. Different projects, done at different times, with different protocols, different machines, different library preps etc.

We have tested combining GTEx and TCGA data for several different purposes where one might expect the analyses' internal controls would have allowed compatibilty, but even then we have always found a persistent batch effect which we were unable to remove.
> 1 votes

# Answer
I am not sure how appropriate it is, but GTEx data was used as the normal controls for TCGA data in Comprehensive analysis of normal adjacent to tumor transcriptomes. It was not just a single panel, but a major part of the paper.

Alternatively, there is RNAseqDB:

> Here, we developped a pipeline that processes and unifies RNA-seq data from different studies. Using the pipeline, we have processed data from the GTEx and TCGA and have successfully corrected for study-specific biases, allowing comparative analysis across studies.
> 1 votes

---
Tags: rna-seq, ngs, cancer
---


# Question
Title: How to set new identities for new clusters based on a gene expression in Seurat?
I have a group of cells that have yellow (EYFP) expression and cells that don't. I want to define the yellow and non-yellow genes in two cell identities.

I clustered the cells and then created Violin Plot for EYFP expression. Then I did FindMarkers between clusters that express high EYFP with low EYFP expression. However, this is not a very accurate way to isolate cells of high and low EYFP expression. I need to use ident to separate EYFP +ve and EYFP -ve cells (but don't know the exact syntax). Then FindMarkers between these two new identities.

I need to subset a certain set of barcodes that express a marker gene (in this case Olfm4).

# Answer
> 1 votes
You have to define your `ident` it has to have the structure of a table. You put the cell names as names and clusters as value (for instance cell 1 = YEPF+ cell 2 = YEFP-):

```
pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

```

Then you can run the `FindAllMarkers` function on this object:

```
FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)

```

---
Tags: seurat
---


title: Highly variable elements?

text:
# Question
Title: Highly variable elements?
Which type of genomic elements in the UCSC browser can be used as an example of the highly labile = highly variable = low conservation regions?

# Answer
> 1 votes
Transposable elements (TEs) are the classic example of rapidly evolving regions in eukaryotic genomes. These are often used as a proxy for neutrally evolving sequences, which are by definition not conserved/constrained. The UCSC genome browser includes various repeat tracks (e.g. RepeatMasker) that includes TEs.

Within protein coding regions, 4-fold degenerate sites are often used as a proxy for neutrally evolving sequences. These sites are where any nucleotide at this position leads to the same amino acid and therefore protein (synonymous mutations).

Potentially even more variable sequences would be those under positive selection, but there is not a single class of element that falls into this category i think.

---
Tags: genomics
---


# Question
Title: Why do I have >10,000 cells in the 10X matrix produced by cellranger?
In my scRNA seq experiment, single-cell libraries were generated using the GemCode Single-Cell Instrument and Single Cell 3′ Library & Gel Bead Kit v2 and Chip Kit (10x Genomics) according to the manufacturer’s protocol. About 8700 cells were added to each channel with a targeted cell recovery estimate of 5000 cells.

However, after processing my raw data with CellRanger software, I get 737,280 cells (barcodes) and about 31,000 genes. I cannot understand why I get so many cells and why this barcode whitelist is not up to 10,000 cells. After filtering this data retaining genes expressed in \>= 5 cells as well as retaining cells with at least 200 detected genes, I get around 2200 cells.

Thanks

# Answer
The raw data from cell ranger contains all of the barcodes detected in the experiment. These raw barcodes are filtered by cell ranger to identify the barcodes that likely represent cells rather than empty droplets / dead cells. You can find the filtered data in the filtered `feature-barcode matrices`. You can also examine the `Run summary HTML` file to see a summary of the cell ranger barcode filtering.

The reason why there are so many raw barcodes is an important optimization by 10X to increase the cell capture rate compared to techniques such as drop-seq. In a 10X experiment the number of gel beads is much higher than the number of cells which leads to a higher probability of a cell and a bead ending up in the same droplet. This also leads to a large number of empty droplets with just a bead and some ambient RNA or a dead cell. These all end up being sequenced and the empty drops tend to have a small number of unique molecules (the majority less than 10 UMIs / barcode).
> 3 votes

---
Tags: scrnaseq, 10x-genomics
---


# Question
Title: How to retrieve logical expressions (KO based) for reactions from KEGG?
The completeness of a module can easily be checked by looking at the `Definition` entry associated with the module. For example, in module M00010, it is given as:

```
Definition  K01647 (K01681,K01682) (K00031,K00030)

```

which can be translated to:

```
K01647 AND (K01681 OR K01682) AND (K00031 OR K00030)

```

If this expression evaluates to `TRUE`, the module is considered complete.

I am wondering whether analogous information exists for a single *reaction*. For example, when looking at reaction R00352, one finds the following information about its Orthology:

But I can't easily find the logical relation that correspond to the KOs for this reaction. It could be:

```
K01648 AND K15230 AND K15231

```

or:

```
K01648 OR K15230 OR K15231

```

or:

```
K01648 OR (K15230 AND K15231)

```

and so on.

In the above example, I know that the correct expression is:

```
K01648 OR (K15230 AND K15231)

```

i.e. one either needs either `K01648` or both of the other two subunits.

Can this information be retrieved from KEGG for each reaction and if so, how?

# Answer
> 1 votes
Some information seems to be available on the module level utilising the API.

Reaction R00352 points to the module M00173. While some information is in the `Definition`, it's not trivial to link it to the particular reaction. This linkage is however easier accessible via API:

---
Tags: kegg, database, data-retrieval
---


title: What's a good ontology for drug names?

text:
# Question
Title: What's a good ontology for drug names?
I have ...

* A database with patient phenotypes in, stored as HPO terms
* Genetic data in whatever format I need

I want ...

* To store drug names in a way that won't make my life difficult
* If it's feasible, to make them so that you *could* draw links between HPO terms, medications, genetic changes, and outcomes. *Note that I'm not actively trying to do that, but I don't want to exclude the possibility of it in the future*

Ideally, I'd like those drug names in some kind of ontology, and even more ideally an ontology that's available on http://bioportal.bioontology.org

So - what's the best drug ontology to use in this case?

**Update**

> What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease? Drugs by disease would be the most useful. The information that would be most useful would be something like

* Atenolol
    + Atenolol 25mg
    + Atenolol 50mg
    + Atenolol 100mg
* Propanolol
    + Propanolol 10mg
    + Propanolol 40mg
    + Propanolol 80mg
    + ...

In short, being able to get a list of dosages for a given drug. It'd be extra-useful if I could also retrieve a set of drugs for a given condition, so in the above example to retrieve a list of beta blockers (Atenolol, Metoprolol, Nebivolol, ...)

> Have you looked at any of the various drug ontologies listed in bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

* MDDB doesn't seem to be structured as above
* Likewise The Drug Ontology
* RxNorm seems to be the best available that I've seen

# Answer
> 3 votes
Since this question's been asked there's been a (what looks like a first useful?) release of Drug Ontology DRON. Perhaps worth checking again?

---
Tags: gene, drugs, ontology, human
---


# Question
Title: How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes?
I have a list of human genes that have been proven to be up-regulated in a disease condition of interest, using microarray analysis. I also have a protein which according to literature is said to be a DNA binding protein. I wish to be able to computationally predict if the protein mentioned binds to any of the promoter sequences of my differentially expressed genes and if yes , then at which position?

Are there any Bioinformatic tools available for this task?

# Answer
> 4 votes
You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you binding sites: genomic intervals where those TFs bind.

https://bioinformatics.stackexchange.com/a/2491/776

You can then map or intersect promoter regions with those TF binding sites with BEDOPS `bedmap` or `bedops`, respectively:

https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html

The use of `bedmap` is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.

# Answer
> 0 votes
If you work with human, I would start by intersecting your promoter coordinates with the ReMap2018 database which is a comprehensive collection of published ChIP-seq datasets. Extract those TFs that fit your research question. Alternatively you can browse NCBI for published datasets of your TF and then intersect those with your annotations.

Also cross-posted: https://www.biostars.org/p/381751/

---
Tags: computation, homer
---


title: How I normalize these two sets of data

text:
# Question
Title: How I normalize these two sets of data
I have average log fold change for a cluster of cells versus another cluster of cells by Seurat like below

```
         avg_logFC
TAGLN   -2.507833541
TXNRD1  1.093672244
MYL9    -1.579165463
ACTB    -1.482386539
ACTA2   -2.379757433
SQSTM1  1.227562767
ACTG2   -2.144757086

```

I have also mean log2 ratio of a set of proteins in treatment versus control like below (I converted protein IDs to gene symbol)

```
mean log2ratio  Gene Symbol
-0.3    CYFIP1
0.4     GNPTG
-0.3    WWC2
-0.5    TRPV4
-0.4    CHRM2
0.3     MFSD11
-0.3    ANXA11
0.3     CENPE
-0.4    RPL22L1
0.3     SAMD1

```

How I can normalize these two data or make them compatible for plotting a heatmap, pathway analysis, etc?

# Answer
There's no need to normalize them, you're not comparing them. Just use them as they are.
> 2 votes

---
Tags: r, scrnaseq, normalization, proteomics, heatmap
---


# Question
Title: Confusion related to effective population size
I am making a demographic model for African American population in which I say 80% of ancestry comes from African population and 20% comes from European population. I am using the time/population time from Gazave et al. (2014).

The authors doe not mention the effective population size (Ne) of African American population after admixture. Is there a publication someone can point me to that mentions the effective population size of African American (AA) or even European American (EA). I know there is not much difference between the two. It's just I need a figure to have my simulation up and running.

**Note** *N*e is an estimate of the people size based genetic diversity and migration and is associated with the Fst statistics

# Answer
> 0 votes
That is a fantastic question and in my opinion a really important question in context. The Ne of Africa will be much bigger than the European population, however I'm less sure about the West African population which gave rise to the Afro-American population. I suspect it is still large because of the "out of Africa" theory.

My advice is to start by taking a complete look at the entire subject, from Tenesa et al (2007). Ne is a really important parameter in any population study and extremely useful. I have feeling this could be complex, but I cannot hypotheses any further because I don't do human genetics (except for immunology).

---
Tags: phylogeny, phylogenetics
---


# Question
Title: Do any publicly available databases detail protein structure and functional domains?
I am interested in finding a database that takes a gene or protein name as input (possibly with the option to specify transcript) and gives information about the protein's functional domains in terms of either specific residue ranges, genomic coordinates, etc. Ideally, the database would include 3D conformational information as well as the functional domains. Does such a database exist? It seems that this information can be inferred for a given protein using various publications and references but I have yet to find a single collection containing these data all in one place.

# Answer
> 6 votes
Some of this information (at least some domains, active sites, etc) is available from UniProt.

If you want to download their whole database, you can search without specifying any terms and then click the Download button.

# Answer
> 3 votes
EnsEMBL also has this. Search for your gene of interest, choose your transcript, go to the page of its protein product(s), and select "Domains & Features" from the right-hand menu (using human p53 as an example):

```
Domain source       Start End  Description                                            Accession InterPro
PANTHER             3     331  -                                                      PTHR11447:SF6-
Pfam                5     29   p53 transactivation domain                             PF08563   IPR013872
PANTHER             3     331  p53 tumour suppressor family                           PTHR11447 IPR002117
PRINTS              116   142  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              158   179  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              213   234  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              236   258  p53 tumour suppressor family                           PR00386   IPR002117
Prosite_patterns    237   249  p53 tumour suppressor family                           PS00348   IPR002117
PRINTS              264   286  p53 tumour suppressor family                           PR00386   IPR002117
Pfam                95    288  p53 DNA-binding domain                                 PF00870   IPR011615
Superfamily         97    287  p53-like transcription factor DNA-binding              SSF49417  IPR008967
Gene3D              94    297  p53/RUNT-type transcription factor DNA-binding domain  2.60.40.720IPR012346

```

# Answer
> 0 votes
In addition to UniProt and EnsEMBL as suggested above you might like:

* InterPro, functional analysis of proteins by classifying them into families and predicting domains and important sites.
* ECOD, a classification of protein domains by evolutionary relationship.
* Pfam, a collection of protein families.
* CATH, protein structural domains.
* Genome3D, which combines various resources to structurally annotate proteins.

---
Tags: database, gene, proteins, protein-structure, rna-structure
---


title: PCA vs tSNE in single cell RNA-seq

text:
# Question
Title: PCA vs tSNE in single cell RNA-seq
**What makes tSNE being the preferred dimensional reduction** for visualization in single cell RNA-seq over PCA?

I am aware that tSNE works better at showing local structures and fails to capture global structures of the data.

But I think I don't fully get the reason of **why is this an advantage**? It offers better resolution? Or better separation of the cells as compared to PCA?

# Answer
tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.

As for global vs. local, we are much more interested to see similarity of cells to a limited number of neighbours indicating a celltype and grouping these close together. This is more important than the distance between such celltypes. (assuming the separation in your tSNE is driven by a biologically meaningful factor such as celltype and not some confounder)

Edit:

Since I just made these images for myself anyway I might as well post them. The first 2 PCs show some separation by celltype, tSNE computed over 100 PCs gets very nice separation.
> 11 votes

---
Tags: scrnaseq, ngs, single-cell
---


# Question
Title: Why does IGV turn a VCF QUAL of "." into -10?
IGV is turning QUAL "." into -10 from this VCF of structural variations. I'm wondering if this means anything or if I should report a bug to the IGV devs.

Here's a row from a VCF:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT SAMPLE

chrY    11292098    DUP00008569 C   <DUP>   .   PASS    IMPRECISE;SVTYPE=DUP;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chrY;END=56834684;PE=3;MAPQ=24;CT=5to3;CIPOS=-265,265;CIEND=-265,265    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV0/1:-5.71458,0,-20.0913:57:PASS:0:90774:0:-1:6:5:0:0

```

You see the QUAL value is "."

IGV reports it as -10

(This call is probably a false positive but it's just an example. There are others)

# Answer
Looks like they are treating the "." as a value of 1. According to the VCF spec the quality is "−10log10 prob". Where are you getting your QUALs from? I have never seen the quals as a "." and IGV probably needs to do some translation into an int but a QUAL of -10 makes no sense.

I would report it and let them explain.
> 3 votes

---
Tags: vcf, variant-calling, igv
---


# Question
Title: Removing PCR duplicates in RNA-seq Analysis
After reading some of the forum posts in Biostar and SeqAnswers I find it very confusing whether to filter out the duplicate reads from aligned files or not. As far I understand it's very difficult to distinguish between highly expressed genes and duplicate reads and we may lose important information during the filtration process.

So, is it really necessary to remove the duplicates in differential expression analysis using RNA-seq data?

# Answer
For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them.

In general, the smaller the amount of RNA input into the library prepartion the worst the duplication. Many protocols for very low input quantities (such as single cell) include random barcodes called UMIs (Unique Molecular Identifiers). These allow PCR duplicates to be distinguished from genuinely independent molecules that just happen to have the sample mapping position.
> 9 votes

# Answer
Generally you should just leave them as is. One does remove/mark duplicates in DNA seq.

For further read check this Nature paper
> 3 votes

# Answer
Some of the researchers in my institute have researched on this and use UMIs (Unique Molecular Identifiers) to eliminate PCR duplicates. They have been able to increase the reproducibility of the rna seq analysis they had been working on. The UMIs capture some lengths of diverse nucleotides from the rna-sequence. I cannot validate how much UMIs can help in deduplication, as even Nature journal says that UMIs specifically work on the 3’ end. But it has been observed that eliminating duplicates doesn’t improve the accuracy of quantification and it could introduce a bias in the data. So do not eliminate the duplicates, it's just a waste of time.
> 0 votes

---
Tags: rna-seq, differential-expression, quality-control
---


# Question
Title: How do i map subcellular localizations to lower the number of classes to predict
For a project, I want to predict the subcellular localization of some proteins. I am a student in software Informatics at the University of applied Science HEIG-VD.

Subcellular localization prediction is a common problem in the bioinformatics, but I have a problem

With the data I found. In this database: of human and mouse proteins, we can see that there are more than 30 different localizations.

At my level, and for the use of my exercise, 30 classes to predict is too much.

Where can I find a mapping to reduce those 30 classes to 4 or 5?

I tried QuickGO, to see the ancestor tree of the different classes, but nothing conclusive. Somebody might have done that, but I lack knowledge on bioinformatic web research.

There are the classes I wanted to map to reduce the number of classes :

* Cytoplasm (GO:0005737)
* Cytoplasmic Vesicles (GO:0016023)
* Endoplasmic Reticulum (GO:0005783)
* Endosomes (GO:0005768)
* Extracellular (GO:0005576)
* Golgi Apparatus (GO:0005794)
* Lipid Particles (GO:0005811)
* Lysosomes (GO:0005764)
* Melanosome (GO:0042470)
* Mitochondria (GO:0005739)
* Nucleus (GO:0005634)
* Peroxisome (GO:0005777)
* Plasma Membrane (GO:0005886)
* Synaptic Vesicles (GO:0008021)
* Cellular Component Unknown (GO:0008372)
* Apical Plasma Membrane (GO:0016324)
* Basolateral Plasma Membrane (GO:0016323)
* Centrosome (GO:0005813)
* Golgi Cis Cisterna (GO:0000137)
* Cytoskeleton (GO:0005856)
* Early Endosomes (GO:0005769)
* ERGIC (GO:0005793)
* Inner Mitochondrial Membrane (GO:0005743)
* Late Endosomes (GO:0005770)
* Medial-Golgi (GO:0005797)
* Outer Mitochondrial Membrane (GO:0005741)
* Secretory Granule (GO:0030141)
* Golgi Trans Cisterna (GO:0000137)
* Golgi Trans Face (GO:0005802)
* Tight Junction (GO:0005923)
* Transport Vesicle (GO:0030133)

# Answer
> 1 votes
One way to do this would be to plot the first 10 GO categories of the the CC sub ontology (CC from cellular component). Or you can just visit the official webpage to see the list of children terms of them.

However there aren't only 4 or 5 categories, so you could arbitrary pick some of those.

# Answer
> 0 votes
This is probably one of those things where as long as you can justify your choice somewhat.

It might be worth seeing what others have done (see loctree3 and PSORT) and seeing how well they match your case, then add in any subdomains you have a particular interest in.

---
Tags: proteins, gene-ontology
---


title: bedtools feature out of bounds

text:
# Question
Title: bedtools feature out of bounds
I'm trying to count the GC in regions leading up to polyA clusters. I found the bedtools command nuc, one of the outputs of which is the GC content.

I took the mouse polyA cluster database from here, and the mm10 primary assembly soft masked from here. I imported the polyA clusters bed file to R. I then got a list of all the polyA clusters in the genes I'm interested in. For positive strands on polyA I took the start position, took 100 away from it and made that the start of the region of interest, then made the start of the polyA signal the end. For things on the negative strand I added 100 to the end of the polyA signal. I then exported the results to a bed file format. The code I used to do this is below.

When I try to run:

bedtools nuc -fi Mus\_musculus.GRCm38.dna\_sm.primary\_assembly.fa -bed PolyA.bed

Every signal returns some variant of: Feature (chr10:21,080,315-21,080,415) beyond the length of chr10 size (0 bp). Skipping.

I can however look at several of these regions in the UCSC genome browser without a problem, at least I don't think I had a problem.

Can anyone point me to what's going wrong or to another way I can solve this.

```
GetRanges <- function(FilteredClusters){
  for (x in 1:nrow(FilteredClusters)){
    if(FilteredClusters$Strand[x] == '+'){
  FilteredClusters$RegionEnd[x] <- FilteredClusters$ClusterStart[x]
  FilteredClusters$RegionStart[x] <- FilteredClusters$RegionEnd[x] - 100
}
else if(FilteredClusters$Strand[x] == '-')
    {
      FilteredClusters$RegionStart[x] <- FilteredClusters$ClusterEnd[x]
      FilteredClusters$RegionEnd[x] <- FilteredClusters$RegionStart[x] + 100
    }
  }
    ClusterCut <- select(FilteredClusters, 'Chr', 'RegionStart', 'RegionEnd')

    return(ClusterCut)
}

```

# Answer
> 3 votes
There is no chromosome named `chr10`, but there is one named `10`. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of `chr` in your BED file and fix things like `chrM` -\> `MT`.

---
Tags: rna-seq, bedtools
---


# Question
Title: Issues performing variant calling with GATK
I am trying to perform variant calling on a BAM file generated through STAR version **STAR\_2.6.0b** for wheat genome using GATK haplotypecaller as follows:

```
gatk HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
02:34:13.450 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
02:34:13.600 INFO  HaplotypeCaller - ------------------------------------------------------------
02:34:13.600 INFO  HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.0.4.0
02:34:13.600 INFO  HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
02:34:13.601 INFO  HaplotypeCaller - Initializing engine
02:34:40.721 INFO  HaplotypeCaller - Done initializing engine
02:34:48.175 INFO  HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
02:34:48.598 INFO  HaplotypeCaller - Shutting down engine
[June 7, 2018 2:34:48 AM PKT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.59 minutes.
Runtime.totalMemory()=10305929216
java.lang.IllegalArgumentException: samples cannot be empty

```

When I searched for this error the most proposed solution was to validate BAM file as suggested here. The ValidateSamFile gave following output for the bam file:

```
gatk ValidateSamFile -I sorted.bam -M SUMMARY
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar ValidateSamFile -I sorted.bam -M SUMMARY
02:53:27.498 INFO  NativeLibraryLoader - Loading libgkl_compression.so   from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 02:53:27 PKT 2018] ValidateSamFile  --INPUT sorted.bam --MODE   SUMMARY  --MAX_OUTPUT 100 --IGNORE_WARNINGS false --VALIDATE_INDEX true --INDEX_VALIDATION_STRINGENCY EXHAUSTIVE --IS_BISULFITE_SEQUENCED false --MAX_OPEN_TEMP_FILES 8000 --SKIP_MATE_VALIDATION false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 02:54:08     SamFileValidator        Validated Read    10,000,000 records.  Elapsed time: 00:00:36s.  Time for last 10,000,000:   36s.  Last read position: TGACv1_scaffold_080127_1DS:28,522
INFO    2018-06-07 02:54:44     SamFileValidator        Validated Read    20,000,000 records.  Elapsed time: 00:01:12s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_602746_7DL:127,507
INFO    2018-06-07 02:55:20     SamFileValidator        Validated Read    30,000,000 records.  Elapsed time: 00:01:48s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_641946_U:17,161
INFO    2018-06-07 02:55:55     SamFileValidator        Validated Read    40,000,000 records.  Elapsed time: 00:02:23s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_471755_6AL:43,097
INFO    2018-06-07 02:56:31     SamFileValidator        Validated Read    50,000,000 records.  Elapsed time: 00:02:59s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_434075_5DL:29,039
INFO    2018-06-07 02:57:07     SamFileValidator        Validated Read    60,000,000 records.  Elapsed time: 00:03:34s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_062491_1DL:40,466
INFO    2018-06-07 02:57:42     SamFileValidator        Validated Read    70,000,000 records.  Elapsed time: 00:04:10s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_558158_7AL:13,072
INFO    2018-06-07 02:58:18     SamFileValidator        Validated Read    80,000,000 records.  Elapsed time: 00:04:46s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_435811_5DL:29,008
INFO    2018-06-07 02:58:53     SamFileValidator        Validated Read    90,000,000 records.  Elapsed time: 00:05:21s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_572030_7AS:9,114

## HISTOGRAM    java.lang.String
Error Type      Count
ERROR:MATE_NOT_FOUND    107771
ERROR:MISSING_READ_GROUP        1
WARNING:MISSING_TAG_NM  96856336
WARNING:RECORD_MISSING_READ_GROUP       96856336

```

In order to remove the errors shown by ValidateSamFile I used FixMateInformation as suggested here. as shown below(sharing part of file):

```
gatk FixMateInformation -I sorted.bam -O mate.bam -MC
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar FixMateInformation -I sorted.bam -O mate.bam -MC
06:13:13.934 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 06:13:19 PKT 2018] FixMateInformation  --INPUT sorted.bam --OUTPUT mate.bam --ADD_MATE_CIGAR true  --ASSUME_SORTED false --IGNORE_MISSING_MATES true --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 06:13:24     FixMateInformation      Sorting input into queryname order.
INFO    2018-06-07 06:23:26     FixMateInformation      Sorting by queryname complete.
INFO    2018-06-07 06:23:26     FixMateInformation      Output will be sorted by coordinate
INFO    2018-06-07 06:23:26     FixMateInformation      Traversing query name sorted records and fixing up mate pair information.
INFO    2018-06-07 06:23:33     FixMateInformation      Processed      1,000,000 records.  Elapsed time: 00:00:06s.  Time for last 1,000,000:    6s.    Last read position: TGACv1_scaffold_361122_4DS:70,952
INFO    2018-06-07 06:23:42     FixMateInformation      Processed     2,000,000 records.  Elapsed time: 00:00:15s.  Time for last 1,000,000:    8s.    Last read position: TGACv1_scaffold_044571_1BL:711
INFO    2018-06-07 06:23:50     FixMateInformation      Processed     3,000,000 records.  Elapsed time: 00:00:23s.  Time for last 1,000,000:    8s.   Last read position: TGACv1_scaffold_433135_5DL:36,647
 .
 .
 .
 .
 [Thu Jun 07 06:40:50 PKT 2018] picard.sam.FixMateInformation done.    Elapsed time: 27.62 minutes.
 Runtime.totalMemory()=8004304896
 Tool returned:
 0

```

After performing FixMateInformation I again validated the BAM file using ValidateSamFile but the errors are still the same as shown above which shows Mate not found, which means error is not fixed.

How these errors can be fixed? So that I can perform variant calling using BAM file?

# Answer
> 1 votes
The primary error isn't the missing mate, but likely the missing read groups. Use Picard `AddOrReplaceReadGroups` to ensure your BAM files have read groups and then run those through the haplotype caller.

# Answer
> 0 votes
You could add read groups:

```
['STAR', '--twopassMode', 'Basic', '--genomeDir', '/data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/',
         '--runThreadN', '10', '--limitSjdbInsertNsj', '5000000', '--outSAMtype', 'BAM', 'SortedByCoordinate', '--twopass1readsN', '-1',
         '--sjdbOverhang', '100', '--readFilesIn', sra1 + '_1.fastq', sra1 + '_2.fastq',
         '--outSAMattrRGline', 'ID:' + rg1, 'SM:' + rg1, 'PL:ILLUMINA'], shell=False)

```

# Answer
> 0 votes
The problem is that you are missing read groups. You can add them directly to bam files you already have using picard tools. Something like

```
java -jar ./software/picard.jar AddOrReplaceReadGroups \
  I=input.bam \
  O=output.bam \
  RGLB=lib1 \
  RGPL=illumina \
  RGPU=unit1 \
  RGSM=hello

```

---
Tags: rna-seq, variant-calling, gatk
---


# Question
Title: processing fastq files using cellranger on linux
I am trying to perform a cellranger count on fastq files generated from a 10x genomics single cell RNA Seq run.

Just to provide some background, I ssh’d into the AWS using our AWS IP with MobaXTerm linux command line.

My files are in the Isilon. The file folders are `path/to/runs/` \- this folder contains the `.fastq` files for each sample (ex: `SI_GA_A1_S1`, `SI_GA_B1_S2`, `SI_GA_C1_S3`, and `SI_GA_D1_S4`). I can see these folders listed in the linux command line.

The cellranger is found when I write the following:

```
path/to/cellranger-3.0.1/cellranger-3.0.1

```

then `ls` brings up 13 items located in the folder, one of them being `cellranger`.

When I write the following command,

```
/path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger count \
    --id=run2\
    --transcriptome=/path/to/refseq/refdata-cellranger-hg19-3.0.0\
    --sample=sample1\
    --expect-cells=5000

```

and it gives me:

```
-bash: /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger: No such file or directory

```

This looks correct from what the 10x website says for single cell RNA Seq. Please advise.

```
$ ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger
-rwxrwx--- 1 1003080 1000000 460 Dec 5 19:41 /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger

```

# Answer
> 1 votes
this has been resolved. I downloaded the windows file and this caused trouble but when downloading the file within the Linux, it seems to be working.

# Answer
> 0 votes
It could be that `cellranger` does not have execute permissions. Try

```
chmod u+x [insert internal isilon path here]/[insert windows username here]/cellranger-3.0.1/cellranger-3.0.1/cellranger

```

(or the Windows equivalent, which seems to be to give Execute permission using the GUI)

---
Tags: linux, bash, cellranger
---


# Question
Title: ipfam or pfam domain-domain interactions
I'd like to download domain-domain interactions that were in ipfam (which doesn't seem to be up-and-running anymore).

for an example of what I would like to be able to download via ftp or another source, see this list of interactions with Death domain: https://pfam.xfam.org/family/PF00531#tabview=tab8

Is this available in one of these files in the current release?? ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current\_release

# Answer
> 0 votes
I think you're looking for ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current\_release/database\_files/pfamA\_interactions.txt.gz. Based on the file's name, that's a list of domain interactions. The file has two columns, so presumably each line represents an interaction.

Here's what I get for the Death domain (PF00531):

```
$ zgrep PF00531 pfamA_interactions.txt.gz
PF00036 PF00531
PF00531 PF14786
PF00531 PF00531
PF00791 PF00531
PF01335 PF00531
PF13499 PF00531
PF13833 PF00531
PF14786 PF00531

```

Those are 8 interactions, which matches what the PFam link you gave in your question shows. If you want to convert the accessions to something more informative, you could try Pfam's REST API and do something like:

```
$ zgrep PF00531 pfamA_interactions.txt.gz |
    while read a b; do
    name1=$(wget -O - "https://pfam.xfam.org/family/${a}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    name2=$(wget -O - "https://pfam.xfam.org/family/${b}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    printf "%s\t%s\t%-20s\t%-25s\n" "$a" "$b" "$name1" "$name2";
done
PF00036 PF00531 EF hand                 Death domain
PF00531 PF14786 Death domain            Tube Death domain
PF00531 PF00531 Death domain            Death domain
PF00791 PF00531 ZU5 domain              Death domain
PF01335 PF00531 Death effector domain   Death domain
PF13499 PF00531 EF-hand domain pair     Death domain
PF13833 PF00531 EF-hand domain pair     Death domain
PF14786 PF00531 Tube Death domain       Death domain

```

Obviously, if you need anything more fancy than that, it's time to write a proper parser using Perl or Python or any other language that can understand XML.

---
Tags: interactions, domains
---


# Question
Title: How do you normalise read coverage in a BAM file?
This is a question from the Oxford Nanopore community, from user *Michael Radzieta*:

> I have sequenced some plasmids using the rapid barcoding kit, I have attempted to assemble the data using several programs (unicycler, spades, ra) but nothing has been able to assemble the plasmids I am expecting. I have read that if the coverage is too high, then it can negatively affect the assembly. Based on the approximate size of my plasmids and the number of filtered reads I have, the coverage is estimated at \>5000.
>
> Can anyone recommend a tool to subsample or normalise nanopore reads to a certain depth, such as 100x coverage?

# Answer
> 2 votes
I just updated one of my SAM-normalisation tools to more accurately normalise mapped reads to a target coverage. It uses reservoir sampling on a SAM-formatted file to choose reads to bring the coverage up to 100X (if possible).

https://gitlab.com/gringer/bioinfscripts/blob/master/samNormalise.pl

Example Usage:

```
samtools view -h input.bam | ~/scripts/samNormalise.pl | samtools view -b > normalised.bam

```

The output BAM file can be converted into a fastq file using 'samtools fastq':

```
samtools fastq normalised.bam > normalised.fastq

```

Here are some example before and after shots \['samtools mpileup input.bam'; 'samtools mpileup normalised.bam'\]:

Notice that where the coverage is below 100, the coverage is identical. Coverage is below 100 in some areas because some reads end, and no additional reads start until later on in the reference sequence.

The options are coded as variables at the beginning of the script at the moment, but I can easily convert them to command-line options in the future if there's demand for it.

This will work for long reads (i.e. PacBio, Nanopore), bearing in mind that there may be more low-coverage gaps if the reads are sparsely placed.

As implemented, this won't work properly with paired reads and will leave behind singletons.

---
Tags: sam, nanopore
---


title: ANNOVAR, 1000g database calling

text:
# Question
Title: ANNOVAR, 1000g database calling
I have recently downloaded Annovar, and succesfully filtered my files with a few databases, downloaded from Annovar. But when I try to do the exact same thing with the 1000 genomes project database, I get this error many times in the terminal, one for each line of \**.

```
Argument "&" isn't numeric in numeric eq (==) at annotate_variation.pl line 1962, <DB> line &&.

```

`&` can be `A`, `C`, `G`, or `T`.

`&&` every line of the database

I have tried filtering with `1000g`, `1000g2010` and `1000g2010jul` databases, but all give these errors.

An example of the commands I use for filtering:

```
perl annotate_variation.pl --filter --dbtype CEU.sites.2010_07 testfile1.txt humandb/

```

Can anyone tell me if I'm doing something wrong, if this is known, and if I can do anything about it?

# Answer
That's a warning. Presumably, a match failed (`$&` is the default Perl variable which holds the result of a match operation). So something in your input is not exactly what the script was expecting. You can probably safely ignore these messages.

Alternatively, update your question and add a few example lines of input which we can use to reproduce the error and maybe we can figure something out. If that isn't possible, at least show us line 1962 of the `annotate_variation.pl` script since that's the one that's failing.
> 1 votes

---
Tags: sequence-annotation
---


# Question
Title: Store and search on variant and annotation data
I need to store a large number of variants with their annotations. I need to import dbSNP, as well as Clinvar, and ExAC etc. Also I have sample variant data from patients. I need to be able to search on genomic regions as well as annotations. And it needs to be scalable. What would be a preferred data storage technology for such a case? What is the industry using? It needs to be accessible by APIs/REST. We are using ElasticSearch, but that is proving to be expensive in terms of servers needed.

Something similar to this has been discussed in The state, limitations and comparisons of large variant stores, but that post is a little old, I am sure there has been progress since then.

# Answer
The question you linked is about storing genotypes. Your question is about keeping site information, which is different and much easier.

As to solutions, you can parse the information you want and store them in a standard SQL database. You can use UCSC binning for region queries. Modern SQL databases all provide R-tree indices, which can be adapted for region query and is probably better these days.
> 1 votes

---
Tags: database, snp, indel
---


# Question
Title: Kruskal–Wallis and wilcox.test choice of test for expression data
I have bit of confusion regarding how ggpubr has implemented the KW and wilcox.test. By default usage of ggpubr reports p value for wilcox.test() but when it comes to Kruskal–Walis it is different .

But as i have to compare more that two groups im opting for Kruskal–Wallis so now the confusion is i get a p value which is significant but if i do the test independently then the p value is not significant not sure why is that I have asked similar question old question

My data

```
dput(df1)
structure(list(mHSC = c(7.7559725, 3.307901197229, 7.91049997141824,
10.9659075, 5.66225942792448, 9.40594075, 8.03212169895329, 10.929870369027,
9.80897606357472, 4.29298594445745, 8.58579242742772, 12.0936123426199,
8.61371460307564, 11.0174749823542, 10.5147718423478, 5.53565725,
8.12612925, 10.6881344317581, 8.48826962785442, 5.52294757621421,
8.38046925, 9.44552809599148, 9.19142493299942, 7.03072998196466,
8.88679975, 0.3126697975, 10.4907798679713, 7.864884, 12.7182883447252,
6.3180225, 5.87640020441158, 7.87333318550091, 8.15221192095875,
4.51086624773858, 5.957332, 6.90057575, 5.64708909092112, 10.7124580115802,
9.82466018308675, 8.49078, 8.63046025, 6.96721675, 7.57871225,
8.815752, 6.72449933154642, 10.7951750189282, 8.01780331997795,
6.07469990935439, 8.31037653431384, 9.60791091202699, 7.11738304581177,
6.27364283988744, 10.8712275, 4.76963356477331, 8.22113922937348,
4.68111125, 7.81847825, 8.4123615, 7.1236546846872, 5.5773485,
10.72958, 7.716703, 7.606539, 9.3070835272045, 7.75540042667672,
9.8687305, 8.98201221147182, 8.93182675, 10.868215, 10.8990875,
10.9377104128057, 8.60048104950787, 7.81762687302046, 9.24377416902964,
6.57057975, 7.52549392459474, 5.9971085, 10.7335671888057), mCMP = c(8.08511525,
3.15081874988432, 8.21803244045432, 12.7012625, 6.95234868015311,
10.24927775, 8.02201688082928, 11.1415993774267, 8.6576364168622,
4.33979841956685, 9.20120340882164, 12.6360469440726, 8.82394480426142,
11.7614369660881, 10.518813583537, 5.9224525, 8.3939755, 10.9519826793289,
9.07562352307607, 4.68003855905254, 8.7914695, 9.35721245249627,
9.52072790562384, 7.52153767036027, 9.37291475, 5.323077, 11.316187036788,
8.08079325, 12.9250246147509, 6.3281185, 6.36628784788811, 8.26354221670866,
8.80757801591239, 5.08139074215397, 6.6258195, 7.56676875, 6.30618317691088,
12.8268812686057, 10.4361513118013, 9.12211125, 8.94618675, 7.25775925,
8.7313635, 9.8282665, 7.10318403743241, 10.9016853993601, 8.75786307748107,
6.2613305781586, 7.24503928286349, 10.1900260653741, 6.04352467875968,
5.41919427988363, 11.689205, 5.2447156735349, 7.81023802636728,
5.484589, 7.9321635, 8.8598865, 7.26208235921707, 11.7873, 11.256575,
8.574748, 8.3324525, 9.69726624009281, 8.51491553090231, 10.792235,
9.97113940283311, 10.08911925, 11.867595, 11.8185525, 11.9422571106972,
9.19311951548278, 8.61818925695215, 9.34701182556221, 7.063838,
7.41275116679947, 6.84360325, 10.478012600725), mGMP = c(8.07654075,
3.12607874540186, 7.91384610943697, 12.9717825, 6.533351883655,
10.394515, 8.26093661667729, 11.2346407524235, 11.1019372952288,
4.34824706553879, 9.19896153987714, 12.6475029547774, 8.53466968298702,
11.6750087661439, 10.7705094222618, 5.60385575, 8.32733375, 10.6929416341222,
8.85393492258581, 4.40609754458633, 9.2379455, 9.43854664655775,
9.66597253836355, 7.99640783151295, 9.60382225, 0.40676415, 11.3480173068737,
8.2298565, 13.3988341481968, 5.2305385, 6.12320931339542, 8.16212407932226,
8.98321413854527, 5.21894298832466, 6.88443175, 7.9935035, 6.66952804861909,
11.3240966560241, 10.795029283416, 8.830129, 8.58251875, 7.6113605,
8.508246, 9.63661975, 6.87344684043907, 11.2237814311451, 9.43948359781618,
6.10350570815157, 8.37769527299503, 9.98272594007505, 7.15195399920282,
5.48458261850125, 11.8302525, 5.66643624946888, 8.36658108147322,
5.9928725, 8.56603575, 9.670526, 7.70714893990735, 9.552538,
11.5112725, 8.99551975, 9.504391, 9.76913665769114, 8.58112436662299,
11.2680775, 10.5535730280208, 10.226635, 11.9211375, 12.0610625,
12.3249064903557, 9.52625419135532, 8.66295412469679, 10.0469451814325,
7.700204, 7.44384030300093, 7.40615025, 10.222430040369), mMono = c(6.86074675,
0.899186661852538, 5.72102113412742, 11.0054, 5.22275722758541,
8.56569375, 7.25543228485158, 9.37557140457706, 7.92766781596895,
1.22333714157354, 7.54589243251616, 11.0881258432345, 7.91303390155615,
10.3482735440588, 7.2731491579751, 3.10091425, 6.91085225, 9.03107467062931,
7.91879733419338, 2.71978845710368, 8.31396425, 6.89137159431571,
8.68813944060236, 6.65073272892769, 8.41960525, 0.283334725,
10.0529313221945, 6.7606855, 11.6800576702458, 3.6173355, 3.81146596960356,
6.93219836413735, 7.92727013246883, 4.09314670396697, 5.78634575,
6.55346475, 5.05937426539404, 10.0514401266737, 9.21573912287031,
7.80863325, 7.80450125, 5.96109, 6.07350675, 7.14264825, 5.20991063841094,
8.9541923959891, 5.48494732937368, 3.46755807112085, 5.91011468781522,
6.19799161521768, 4.09986325310833, 1.58331376348435, 11.021745,
4.62557956951436, 7.48457597125806, 4.9819695, 7.567706, 8.547291,
7.13353425765299, 9.4336955, 10.6991475, 7.3555725, 8.32300975,
9.14731426222114, 8.2750379299376, 10.85821, 9.63317457338075,
9.6021585, 11.4682025, 11.4990425, 11.4153866521691, 9.02074441245752,
7.96069342402473, 9.45598020603536, 6.462563, 6.78696488936524,
5.50898575, 8.79251593030914)), class = "data.frame", row.names = c(NA,
-78L))

```

Code for the first figure with wilcox test and KW as argument

```
  head(df1)

ex <- melt(df1, id.vars=c(NULL))

head(ex)
my_comparisons <- list( c("mHSC","mCMP"),c("mCMP","mGMP"),c("mGMP","mMono"),c("mHSC","mMono"),c("mHSC","mGMP"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons,method = "kruskal.test")+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

```

New figure

The one im not able to get is in the first figure even with Wilcox and KW as argument the global p value for KW doesn't change where as the KW p value for comparison gives p value which is not significant .

Am i doing it correctly ?

# Answer
The lower graph doesn't look good. The "Wilcoxan values" are all identical and I can't see how if all pairwise comparisons were 0.46 that would result in a KW of \<\<0.001. The way the values hold together in the upper graph looks fine.

Just to summarise

* KW = 1 way ANOVA ... non-parametric version
* Wilcoxan = T test ... non-parametric version

Looking at your code you specify the pairwise comparisons but then in the ggplot2 code you state `method = "kruskal.test"` ... are you sure that isn't `method = "wilcoxan.test"`?

KW will want to perform a global comparison across multiple single columns of data and this might have resulted in some sort of weird mean value being generated, because it will just want to output a single value (my bet is that is why all values are identical). Wilcoxan will be expecting to a request to perform pairwise tests.

---

Okay my advice is to perform an independent Wilcoxan test outside of ggpubr

```
wtest <- wilcox.test(mCMP ~ mHSC, data = df1, paired = TRUE)
wtest

```

**Output** Wilcoxon signed rank test data: weight by group .....

```
wtest$p.value

```

**Output** p value

I would then simply remove the calculation from with the boxplot and present the p values independently. *The only other possibility is whether the data is paired or not and you may have changed that setting between calculations*.
> 1 votes

---
Tags: r
---


# Question
Title: Enrich for genes expressed in the plasma membrane
I have a list of cell markers differentially expressed between 2 groups of cells, that I got with the FindMarkers command in Seurat. I want to enrich for those genes that are expressed in the cell membrane. Are there any specific tools/commands I can use to do that? Thanks!

# Answer
You can do functional enrichment analysis, for instance with `gProfiler` in R and see if your marker genes are enriched in any GO term related with the plasma membrane.

Alternatively, you can generate a list with genes whose product are located at the plasma membrane and perform Gene Set Enrichment Analysis, for instance with `fgsea` in R.
> 1 votes

---
Tags: seurat
---


# Question
Title: How does MinKNOW classify 1D reads as "pass" or "fail"?
I have found a couple of sources<sup>,</sup> that indicate that a read in a 1D² run is classified by MinKNOW as "pass" and put into the `fastq_pass` folder if both of the following conditions are met:

1. Both strands were read
2. The average Phred score of all bases in the read was \>9

(Thus a read in a 1D² run goes into the `fastq_fail` folder if either its Phred is too low *or* if the second strand didn't follow the first through the pore.)

However, I have been unable to track down any source indicating how the pass/fail classification works for 1D runs (where only one strand is supposed to be read per read). One might reasonably guess that it's based solely upon whether the mean Phred score is above 9, but inspection of the data from a 1D run reveals that this is not the case. While all "pass" reads seem to have a mean Phred score of at least 9 (the lowest I see is 9.45), some fail reads have scores above 9 too (with the highest one being 12.9):

Thus it looks to me like this is not the condition used to classify 1D reads. But then what is?

# Answer
> 6 votes
The classification is indeed based purely on comparing a mean to a threshold... but it's *not* the mean of the Phred scores that represent base quality on a logarithmic scale, but rather the mean of the underlying quality values prior to being transformed into Phreds.

If I compute the "mean" quality for a read in this way - by converting each Phred in the read back to a probability, taking the mean of those probabilities, and then converting that mean probability to a Phred - I find that whether a read is a pass or fail is determined by comparing that score to a threshold of 7:

A mean quality score computed in this way is also included in the `mean_qscore_template` column of the sequencing summary file. If you open that file in a spreadsheet program (it's just a TSV) and sort by `mean_qscore_template`, you will find that the `passes_filter` column is always `FALSE` for `mean_qscore_template` values below 7 and always `TRUE` for values above 7.

(Hat tip to Wouter De Coster for successfully guessing that this was the answer in the comments and pointing me on the right path.)

---
Tags: nanopore, minion, minknow, base-calling
---


# Question
Title: Global or local alignment on sequences for which it is assumed that they have common ancestry
As the title says I would like to know which alignment method(global or local) should I use for sequences for which is assumed that they have common ancestry?

The sequences are 21 aa long peptide chains (proteins).

About the background of the sequences, it's only stated that it's ASSUMED that they are from the same ancestry.

The thing that we assume that they are from the same ancestry makes me confused. I know that on similar sequences with the same length usually, we use global.

# Answer
My initial reaction is that 21 amino acids is too little information to demonstrate these proteins are homologues (which is the question behind your question). A homologous site means it shares a common ancestry and is distinguishable from random chance or convergent evolution. A common ancestor can be demonstrated using phylogeny (trees).

If a regular alignment program cannot align these sequences, such as clustalo or muscle, the assumption they are homologues cannot be easily tested using phylogeny. That would then become complicated at which point you would need to ask whether it is worth pursing the question.

My recommendation is to perform a simple protein-protein blast and assess the e-scores to get some feel about you data. Your question is achievable, but you would solve this in context to additional sequence data. Whether that sequence data is a denser sample or else contiguous sequence data needs further exploration.
> 1 votes

# Answer
If you are trying to figure out the proteins' evolutionary history (e.g. if they are homologous) then you want a global alignment. In fact, there are very few cases where a local alignment is better, local is usually done only because it is far more efficient and less computationally expensive than global. But with tiny sequences of just as few residues such as the ones you have, you would want global almost every time.

The only case I can think of where you might prefer local alignment is when you have sequences that share little overall sequence similarity but you suspect they carry some sort of shared motif. For example, in proteins, that could be a protein domain. In that case, you would want to do a local alignment (possibly using something like blast) to identify those sub-regions of higher similarity.

Note that simple sequence similarity isn't really enough for you to decide whether they are indeed homologous. A more sophisticated approach such as psi-blast or HMMER will give better results. Nevertheless, a simple alignment will help you at least discard some obviously non-homologous sequences.
> 1 votes

---
Tags: sequence-alignment
---


title: Find overlap between VCF files

text:
# Question
Title: Find overlap between VCF files
I have two VCF files and I want to compare the missing rate in each of them, but I want to only look at sites that are present in both files. How could I go about getting a list of positions that are included in both VCF files and then filtering those VCF files to only include those positions?

# Answer
Something like vcf-isec? Quoting:

> Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.
> 2 votes

# Answer
The bedtools intersect command is my go-to for operations like this. From the manual:

> By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. `bedtools intersect` allows one to screen for overlaps between two sets of genomic features. Moreover, it allows one to have fine control as to how the intersections are reported. `bedtools intersect` works with both BED/GFF/VCF and BAM files as input.
> 3 votes

# Answer
You can try something like this in `bcftools`, which is a bit quicker than VCFtools, as it uses block gzipping and indexing.

```
# make sure both files are indexed and in bgzip format
bcftools index vcfA.vcf.gz
bcftools index vcfB.vcf.gz

# find overlapping positions between the two files
bcftools isec -n~11 -c all vcfA.vcf.gz vcfB.vcf.gz > sharedPositions.txt

# filter vcf files for the overlapping positions
bcftools view -T sharedPositions.txt vcfA.vcf.gz -Oz > vcfA.sharedPositions.vcf.gz

bcftools view -T sharedPositions.txt vcfB.vcf.gz -Oz > vcfB.sharedPositions.vcf.gz

```

The perhaps use something like plink2 to look at the missing rate per variant and individual for each file:

```
plink2 --vcf vcfA.vcf.gz --missing --out vcfA.missingStats

```
> 2 votes

---
Tags: vcf, variant-calling, variation
---


# Question
Title: How can I predict if the binding site of a DNA-binding protein is located on the promoter sequence of a gene?
Is there any computational tool or any other method to predict if the binding site of DNA binding protein is actually located on a promoter sequence?

# Answer
You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you predicted binding sites: genomic intervals where those TFs bind, with some statistical significance.

https://bioinformatics.stackexchange.com/a/2491/776

You can then map or intersect promoter regions with those TF binding sites with BEDOPS `bedmap` or `bedops`, respectively:

https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html

You can define your promoter regions however you like, e.g., from the TSS to 2.5kb upstream of the TSS. It's up to you and your experiment.

The use of `bedmap` is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.
> 1 votes

---
Tags: bioconductor, genome, proteins, computation, transposable-elements
---


title: Quality GL entries in vcf file

text:
# Question
Title: Quality GL entries in vcf file
I have a `vcf v4.2` file with structural variants, the 8 columns describing the variant are clear, but in the genotype information I got confused by values of `GL` entries. In the vcf specification they mention that `GL` are likelihoods of the individual genotypes, but as far as I am concerned, likelihood should be a number between 0 and 1, therefore I am suprised to see negative numbers. Then I thought that it might be log-likelihood, but then it's strange that the higher the number is the better. Here is a sample of 5 variants with various called genotypes

```
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-22.4949,0,-4.49485:45:PASS:20:1:11:0:1:4:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-3.1888,0,-43.9892:32:PASS:16:26:13:2:9:2:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-0.190635,0,-49.2888:4:LowQual:26:51:31:2:10:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/0:0,-1.01435,-75.899:10:LowQual:18:54:26:2:14:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   1/1:-19.198,-1.5031,0:15:PASS:70:0:58:0:0:5:0:0

```

The `vcf` file was generated by delly genotyper. I followed their recommendations for germline SV calling. However, negative likelihoods are even in the `vcf` specification, therefore it seems that it's not something `delly` specific.

So, what exactly is `GL` in `vcf` files? And how do they relate to other entries `GP` (genotype posterior probabilities, not in the file by delly) and `GQ` (conditional genotype quality)?

# Answer
> 4 votes
The VCF spec clearly states that these are indeed log likelihood values:

> • GL : genotype likelihoods comprised of comma separated floating point log10-scaled likelihoods for all possible genotypes given the set of alleles defined in the REF and ALT fields. In presence of the GT field the same ploidy is expected and the canonical order is used; without GT field, diploidy is assumed. If A is the allele in REF and B,C,... are the alleles as ordered in ALT, the ordering of genotypes for the likelihoods is given by: F(j/k) = (k\*(k+1)/2)+j. In other words, for biallelic sites the ordering is: AA,AB,BB; for triallelic sites the ordering is: AA,AB,BB,AC,BC,CC, etc. For example: GT:GL 0/1:-323.03,-99.29,-802.53 (Floats)

With log-likelihood values, the higher the number, the better the value. So -2 is better than -3.

---
Tags: vcf, structural-variation
---


# Question
Title: how can i download rna-seq data sample for starting rna-seq analysis by R?
I am looking for an induction into RNA-seq analysis using R. I have yet to decide a data set, or a subject area. Any help in this regard would be welcome.

# Answer
You must have found the airway dataset by googling "RNASeq R tutorial"

https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html

Why won't this work for you?
> 2 votes

---
Tags: r, rna-seq
---


# Question
Title: Are duplicate variants against the VCF standard?
I have been trying to understand an error that the EBI's vcf\_validator gives when run on my vcf file. Consider this minimal file:

```
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   0.00    .   .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

```

If I then run vcf\_validator on it, it complains about:

```
$ cat $(LC_ALL="C" ~/vcf_validator_linux  -i test.vcf |& awk '/written to/{print $NF}' )
According to the VCF specification, the input file is not valid
Error: Duplicated variant chr12:529516:C>T found. This occurs 2 time(s), first time in line 6.

```

This is because the two variants are indeed the same. Changing `AACAC` to `AATAC` at position 529514 just means change `C` to `T` at position 529516.

OK, but how is that against the VCF standard? What makes that an invalid VCF file? I am trying to automatically validate input VCF files and stop the analysis if they're invalid, but this seems like a perfectly valid file. Why is `vcf_validator` not letting it through? The `vcf-validator` perl script from vcftools doesn't complain about this. So is it just a bug in EBI's `vcf_validator` or is it really somehow illegal to have duplicate variants in a VCF file according to the specs?

# Answer
No. The only point the specification makes is that:

> It is permitted to have multiple records with the same POS

I suspect this is *intended* to let the file list different variants occurring at the same position, but this incidentally allows duplicate variants.

This said, you probably don't *want* duplicate variants, because it will mess with any statistics you apply over the variants, and will unnecessarily clutter the file and slow down processing, if this happens a lot. You might also get variant annotations on one variant but not another, which can become a problem depending on how you query your variants.
> 7 votes

---
Tags: vcf, variation, vcftools
---


# Question
Title: Finding the optimal resolution parameter in Seurat's FindClusters function
In `Seurats`' documentation for `FindClusters()` function it is written that for around `3000` cells the `resolution` parameter should be from `0.6` and up to `1.2`. I am wondering then what should I use if I have `60 000` cells? How to determine that?

# Answer
Assuming you have an informative selection of variable genes from which you have constructed a number of **useful** PCs, I'd run a number of iterations with `FindClusters()` as described in the other answer, then choose a level which **overclusters** the dataset (for example, clusters that are visibly separate on a t-SNE or other dimensionality reduction plot should definitely have their own number):

```
seuratobject <- SetAllIdent(seuratobject, id='chosen.resolution')

```

Then run:

```
Seurat::BuildClusterTree()
Seurat::FindAllMarkersNode()

```

Assessing the cluster markers for each node will hopefully give you a good idea on which clusters should be combined. Then you can "combine" the clusters and re-label the cells using something like:

```
library(plyr)
cell.labels <- seuratobject@ident
cell.labels <- mapvalues(cell.labels,
  from=0:16, # cluster numbers
  to=c('A', 'B', 'C', 'C', 'D', 'E', 'E', ... )) # etc

seuratobject <- AddMetaData(seuratobject, cell.labels, 'Combined.clusters')

```

The usefulness of the clustering will very much depend on the selection of variable genes, therefore, depending on the (diversity of the) dataset, you will want to **experiment with selection parameters** or **subset the dataset** and repeat the above procedure.
> 7 votes

# Answer
That is a very general recommendation. Depending on your experiment, you can get a very different number of clusters with the same number of cells at the same resolution.

You can actually use a vector of different resolutions and see which one performs best:

```
pbmc_small <- FindClusters(
  object = pbmc_small,
  reduction.type = "pca",
  resolution = c(0.4, 0.8, 1.2),
  dims.use = 1:10,
  save.SNN = TRUE
)

```

There is no software tool that will tell you what is the best number of clusters. You will have to check the expression of known genes or cluster markers to determine which clusters make the most biological sense. This step is probably the most difficult part of single-cell analysis.
> 3 votes

# Answer
Have a look into clustree to assess the different clusters by clustering them and see different levels.

### Example:

```
library(clustree)
data("iris_clusts")

# plot the tree of clusters
# K1 to K5 correspond to different clustering resolutions
clustree(iris_clusts, prefix = "K")

```

### Using the stability index to assess clusters

The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions. The stability index is automatically calculated when a clustering tree is built. Note that each level of clustering corresponds to a different resolution.

```
clustree(iris_clusts, prefix = "K", node_colour = "sc3_stability")

```
> 3 votes

---
Tags: r, single-cell, seurat
---


# Question
Title: Run Nextflow with file dependencies inside a Docker container?
I’m trying to run a Nextflow workflow in a custom Docker container. Without Docker, the workflow succeeds. But running it inside the container leads to an error because a dependent file cannot be found.

Here’s a minimal script to illustrate the issue (similar to the basic Nextflow script example):

```
params.reference = "/data/reference/GRCh38.fa"

reference = file(params.reference)
fastq_channel = Channel.fromFilePairs('/data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "$reference" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t "${reference}.bwt"
    """
}

```

Running this inside a container will cause the path of the query, `/some/fastq/`, to be mounted inside the container. However, the path to the database, `/data/reference`, will *not* be mounted. Consequently, I get the following execution error:

> stat: cannot stat 'GRCh38-chr21.fa.bwt': No such file or directory

Without the container, the command succeeds.

---

I understand *why* the error occurs: Nextflow doesn’t track the dependency on `$reference` inside the command, and thus doesn’t mount its path in the container.

My question is, **how do I fix this?** The documentation is silent on this regard, and none of the workflows I can find online address the problem. For example, the GATK4 IARC workflows would presumably have the same issue, and aren’t configured for use with Docker.

Just adding the `reference` as a process input is insufficient, since BWA requires various other index files. Am I supposed to add these all manually as inputs?

# Answer
> 3 votes
The following workaround seems to work but ideally I don’t want to have to change the `script` portion of my process so that it works with Docker.

```
params.ref = 'data/reference/GRCh38.fa'

ref = file(params.ref)
ref_index_dir = ref.parent
fastq_channel = Channel.fromFilePairs('data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel
    file ref
    file ref_index_dir

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "reference/$ref" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t reference/${ref}.bwt && touch ${sample}.sam
    """
}

```

Here, `reference` inside the script section refers to the parent directory name of the reference, and needs to be adjusted accordingly. Having to refer to `reference/$ref` is annoying but at least it works — both inside a Docker container and outside of it.

# Answer
> 1 votes
Process input files need to be declared properly to allow nextflow to create the required docker mounts accordingly. Therefore the `reference` file should be declared in the process input block, eh:

```
input:
 set sample, file(fastq) from fastq_channel
 file reference

```

---
Tags: docker, nextflow
---


title: Regular expression struggle

text:
# Question
Title: Regular expression struggle
These are my files

```
SO_8543_IP_3_ablt_R1.fq.gz
SO_8543_IP_4_ablt_R1.fq.gz
SO_8543_INPUT_1_ablt_R1.fq.gz  SO_8543_IP_5_ablt_R1.fq.gz
SO_8543_INPUT_2_ablt_R1.fq.gz  SO_8543_IP_6_ablt_R1.fq.gz
SO_8543_INPUT_3_ablt_R1.fq.gz  SO_8543_IP_7_ablt_R1.fq.gz
SO_8543_INPUT_4_ablt_R1.fq.gz  SO_8543_IP_8_ablt_R1.fq.gz
SO_8543_IP_1_ablt_R1.fq.gz
SO_8543_IP_2_ablt_R1.fq.gz

```

Im running the below line

```
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38  -1 ${i}_R1.fq.gz | samtools view -bS - >${i%}.bam;

done

```

So if i use run this

```
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

```

I get files as such

```
SO_8543_INPUT_1_ablt
SO_8543_INPUT_2_ablt
SO_8543_INPUT_3_ablt
SO_8543_INPUT_4_ablt
SO_8543_IP_1_ablt
SO_8543_IP_2_ablt
SO_8543_IP_3_ablt
SO_8543_IP_4_ablt
SO_8543_IP_5_ablt
SO_8543_IP_6_ablt
SO_8543_IP_7_ablt
SO_8543_IP_8_ablt

```

The part of the file name that is consistent is `_R1.fq.gz` which goes after the input argument. As these are single end files so only 1 im giving one end as input.

Error message

```
***
Error: Must specify at least one read input with -U/-1/-2
(ERR): bowtie2-align exited with value 1
./Chip_SEQ.sh: line 25: syntax error near unexpected token `done'
./Chip_SEQ.sh: line 25: `done'

```

it looks like its not able to recognize as input .Any help or suggestion would be highly appreciated

# Answer
> 6 votes
You have an extra `done`, instead you want:

```
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}_R1.fq.gz | samtools view -bS - >${i}.bam
done

```

As an aside, I'm not a big fan of piping `ls` to cut and then uniq. Why not instead:

```
for f in *.fq.gz; do
    sampleName=${f%%_R1.fq.gz}
    echo $sampleName
    bowtie2 -q ...
done

```

That will do the same thing, but (A) it's clearer what's going on and (B) it can be reused if your sample names change. Also, stop writing bash scripts and use snakemake or something like that instead, it scales much better and is less error-prone.

---
Tags: shell
---


# Question
Title: Error: sort order of reads in BAMs must be the same
I am a novice to bioinformatics trying to run a cuffdiff job in cufflinks and keep getting this error message:

```
Error: sort order of reads in BAMs must be the same

```

**GENERATING SAM FILES:**

```
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_B73_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_Zea_mays.AGPv4.dna.toplevel -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_2.fq

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11R-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_Mo17_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2//Mol17_Ref_Lai -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_2.fq

```

**CONVERTING SAM TO BAM:**

```
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load samtools
samtools view -bS a_B73_11-5_accepted_hits.sam > a_B73_11-5_accepted_hits.bam
samtools sort a_B73_11-5_accepted_hits.bam -o a_B73_11-5_accepted_hits.sorted.bam

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5
module load samtools
samtools view -bS a_Mo17_11-5_accepted_hits.sam > a_Mo17_11-5_accepted_hits.bam
samtools sort a_Mo17_11-5_accepted_hits.bam -o a_Mo17_11-5_accepted_hits.sorted.bam

```

**CUFFDIFF JOB**

```
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load cufflinks/2.2.1
cuffdiff -b /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5 -u -N -p 8 /fs/project/PMIU0157/Claire_Stoll/Zea_mays.AGPv4.36.chr.gff3 /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5/a_B73_11-5_accepted_hits.sorted.bam /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5/a_Mo17_11-5_accepted_hits.sorted.bam

```

# Answer
I think cuffdiff is not quite giving you the right error...it looks like the crux of your problem is you aligned each sample to a different reference.
> 3 votes

# Answer
As @swbarnes2 noted, the error message is because the chromosome names, sizes, and orders differ between the two reference genomes/species you aligned against. Even if they were the same, you're comparing one sample against another between organisms, which is innately not something likely to work well. Further, I wouldn't expect the Zea Mays annotation to be very useful with whatever `Ref_Lai` is.
> 2 votes

---
Tags: sequence-alignment, sam
---


title: Detecting broad peaks in sRNA-seq data

text:
# Question
Title: Detecting broad peaks in sRNA-seq data
What kind of tool would be appropriate do detect "broad peaks" in small RNA-seq sequencing data?

MACS2 appears to be developed for ChIP-seq data, but I see that there is a `--nomodel` option. Would that make this program usable in my case?

I suppose I should also set the `--shift` and `--extsize` options, but I'm not sure how:

```
  --nomodel             Whether or not to build the shifting model. If True,
                        MACS will not build model. by default it means
                        shifting size = 100, try to set extsize to change it.
                        DEFAULT: False
  --shift SHIFT         (NOT the legacy --shiftsize option!) The arbitrary
                        shift in bp. Use discretion while setting it other
                        than default value. When NOMODEL is set, MACS will use
                        this value to move cutting ends (5') towards 5'->3'
                        direction then apply EXTSIZE to extend them to
                        fragments. When this value is negative, ends will be
                        moved toward 3'->5' direction. Recommended to keep it
                        as default 0 for ChIP-Seq datasets, or -1 * half of
                        EXTSIZE together with EXTSIZE option for detecting
                        enriched cutting loci such as certain DNAseI-Seq
                        datasets. Note, you can't set values other than 0 if
                        format is BAMPE for paired-end data. DEFAULT: 0.
  --extsize EXTSIZE     The arbitrary extension size in bp. When nomodel is
                        true, MACS will use this value as fragment size to
                        extend each read towards 3' end, then pile them up.
                        It's exactly twice the number of obsolete SHIFTSIZE.
                        In previous language, each read is moved 5'->3'
                        direction to middle of fragment by 1/2 d, then
                        extended to both direction with 1/2 d. This is
                        equivalent to say each read is extended towards 5'->3'
                        into a d size fragment. DEFAULT: 200. EXTSIZE and
                        SHIFT can be combined when necessary. Check SHIFT
                        option

```

In my case the tags are the same things as the fragments, so this "extension" thing is probably not appropriate.

Basically, I would like to call peaks based on the plain real coverage of my small reads. Are there other tools that would be more appropriate?

# Answer
> 3 votes
The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

* ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The `de-novo cluster finding` option as copied from the manual might do the trick:

```
DE-NOVO CLUSTER FINDING
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.

```
* segmentSeq is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.

---
Tags: software-recommendation, hts, peak-calling
---


# Question
Title: SnpEff gene annotation with @ such as PCDHB@
I'm using snpEff to annotate some VCFs and getting some spurious annotations that don't have transcript ids. Some of these annotations have a gene name followed by ampersand such as PCDHB@ and I'm not finding any explanation in the docs.

Any idea what this refers to? Entire gene family?

Thanks.

# Answer
The `@` at the end of a name denotes a "gene cluster in chromosome region". This is mentioned in the official HUGO nomenclature guidelines in the "Characters reserved for specific usage" section:

It's also mentioned in the GeneCards documentation:
> 1 votes

---
Tags: vcf, sequence-annotation, variant-calling, functional-annotation
---


# Question
Title: RefSeq transcript naming - transcript vs gene
I have a dataset composed on numerous RefSeq transcripts with names of the form NM\_###### - essentially standard annotated mRNA accession numbers. I notice that numerous transcripts which have different sequences (different start and stop sites in a gene) all have the same NM\_#### accession number. The interpretation is that they are all from the same gene, but this makes no sense to me. Why should different transcripts with different sequences have the same name? I understand this identifies their gene of origin, but they are still different transcripts.

Here is an example from the .bed file I'm using. Note they have different start and stop locations and even have different strand notations, but all have the same name.

```
chr     start        stop       name            score   strand
chrX    3441730     3443690     NM_001085545    0       -
chrX    3700232     3702192     NM_001085545    0       +
chrX    3750925     3752885     NM_001085545    0       -
chrX    31336368    31338328    NM_001085545    0       +

```

# Answer
No, this is the same sequence. A single RefSeq accession will correspond to one sequence. However, nobody said that sequence needs to be unique in the genome. The specific sequence you are talking about, NM\_001085545.1, is a mouse transcript which can be aligned to many different regions of the mouse's chromosome X.

This is your sequence:

```
>NM_001085545.1
agagaagtctaggcaggactggcagactacagcaatgcctcaaccaagaagcctctgcga
gagcttgttctaggaattttcatcaatggggcttttagtcagcagggtcttgagatgcag
ggattccagtctgttagagccacagccagaagccatagccggagccagctacattcctgg
cagtcgcaagcgaaaaagaaacagtttggaggagttggcaacaagttctaatgttcatgg
acctcaaaaccagggaatgtatccacatcaagttctcaactacatctactggaaaagggt
taagatctcatctaatgatgcttatcaaaacttatttttggacgggcatgatagcgacat
taaaatccgtgctctgggaaaaacatggtgtttacacaaagtatttttatgtcagtcagg
ctactttgctaacattctcaaaggtacttggagagaatcacaccatggtgttataaatct
gatcattaagaatgaggatattgatacccgatctctgcattttgtgtttggtgctttgta
cacggatgcggatttgtcaataacacctctggaagttcctcaagttttggcagcagcatg
cctgcttcgggtggatcaagtaattcagcagtgtgaaggaatcatgaaagaaactatcaa
caggaacactgtgtgctcctattatttggcagcagaaacctatagattaaaagctgtaaa
gacgagatgctttgaatggcttctttgcaatttgatggtacatccaagtgtggcacttta
caaggaagtagacttgaagttgatgtatcttctagcactgtcttctgacttactagtcat
gcaaaaggagattgatgtatataccacactaaaaatatggatgttcctttatcttaatcc
atgctggaacggaaccatgaaacagcttttacaacacgcaaacaactggctttccaccca
catggcatatgttgataacatcagttttcttgaaagtgaagaaggactaatatttcaacc
agtgtttaaaaagctgagatttcagcacattatctgtgacttgacttccacaactattct
tgaacaagatcgactaatacctatggcatggttgtcacccatttacaaacaacagtggtt
gactttgctgcgaacacaagaatatggggtaattggaccacaagttatcaatgaacaaga
acttgaagaatgcaccatgaggtgtggtacaatgatccccaaggatggaagatatacttg
gaagtggtcagttggacgacttggctttcctttacgtgtgacctttaccaggcagtgtgt
aattttaaggcaacggtgtcagaggtgtgatggttctgcttgccacaaccatatccgaaa
tgtcatattcagaataactttggtgtgttttgattccaacaaaagagtaactttcagaaa
gacaacaggttataaaatcctcacctttgaatataaggaggagcaaattgtaatgaaatt
ggatagtgatgttctaaccttccctatgtgtatattctgcaatttcctttttgtaaacct
aggaaatgcagaaaacaagtaatcttatcaattattagactgtaagtattttacaagtat
tgaatggttcatttaatatgaaggctgtgtaatgactacagtgaagaagaccaaaaacaa
attcaatattctttggagtctagtaccaccagaagtgactgcttccctccatcaccttat
ttccatgtatttctgaagaaaacctttacttccaacaaagtatttgtcaaatcaaataaa
cctaatcttctgatttaaatttaccatcgagatccaagcttggaatggaaataagagact
atatttgcaaggtaaaacatctcaatcaatttgttaagtagatttggccctagctaaaat
tgagacaaaataaaatcatgccattgttgaaactgttgtg

```

If you BLAT it against the mouse genome using the UCSC BLAT tool, you will see many hits, most of them almost perfect matches and some of them overlapping the regions given in your bed file:

Your bed file matches the following hits in the image above:

```
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3441731   3443690   1960
NM_001085545.1  1960     1  1960  1960   100.0%  chrX   +     3700233   3702192   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3750926   3752885   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   +    31336369  31338328   1960

```

I don't know why your bed file has all of these various matches for NM\_001085545.1, since you haven't told us what the bed file is supposed to be doing. A bed file is just a collection of regions. For all I know, that bed file is just "a few hits for sequence NM\_001085545.1 in the genome". There is no reason to assume that each region in a bed file corresponds to a gene.

So there's nothing wrong with your file, you're just dealing with a non-unique sequence from the mouse genome. And you don't have multiple sequences with the same accession, you have one sequence, one accession, but multiple hits for that sequence in the genome.
> 4 votes

---
Tags: refseq
---


# Question
Title: UniProt flat file. Interpreting format of SUBCELLULAR LOCATION. What does flag mean?
In creating a parser for the UniProt flat file `uniprot_sprot.dat` (directory) in the manual for the section `SUBCELLULAR LOCATION` is listed

> The format of SUBCELLULAR LOCATION is:
>
> CC -!- SUBCELLULAR LOCATION:(( Molecule:)?( Location.)+)?( Note=Free\_text( `Flag`)?.)?
>
> Where:
>
> * Molecule: Isoform, chain or peptide name
> * Location = Subcellular\_location( `Flag`)?(; Topology( `Flag`)?)?(; Orientation( `Flag`)?)?
>     + Subcellular\_location: SL-line of subcell.txt ID-record
>     + Topology: SL-line of subcell.txt IT-record
>     + Orientation: SL-line of subcell.txt IO-record

What does `Flag` mean?

My best guess at present is that these are the text surrounded by parenthesis, e.g.

```
 Note=Not associated with microtubules (PubMed:21558460).

```

but that does not help in understanding why it is named flag.

While I can use my parser to verify this is correct when this section is working, it would be nicer to have a better understanding of what the word flag means in this context.

I also checked the XML schema but no help there.

---

Per request.

More examples from file.

These where gathered using grep.

```
grep -A 3 -n 'CC   -!- SUBCELLULAR LOCATION'  /mnt/.../uniprot_sprot.dat

```

```
71:CC   -!- SUBCELLULAR LOCATION: Host membrane {ECO:0000305}; Single-pass
72-CC       membrane protein {ECO:0000305}.

2163:CC   -!- SUBCELLULAR LOCATION: Membrane {ECO:0000305}; Single-pass membrane
2164-CC       protein {ECO:0000305}.

6003:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

7039:CC   -!- SUBCELLULAR LOCATION: Cell membrane {ECO:0000305}; Lipid-anchor,
7040-CC       GPI-anchor {ECO:0000305}. Note=In microneme/rhoptry complexes.
7041-CC       {ECO:0000250}.

7409:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

8677:CC   -!- SUBCELLULAR LOCATION: Virion {ECO:0000269|PubMed:9499098}. Host
8678-CC       endoplasmic reticulum-Golgi intermediate compartment
8679-CC       {ECO:0000269|PubMed:15016891, ECO:0000269|PubMed:9499098}.

12871:CC   -!- SUBCELLULAR LOCATION: Nucleus {ECO:0000250|UniProtKB:P48349}.
12872-CC       Cytoplasm {ECO:0000250|UniProtKB:P48349}. Note=Translocates from
12873-CC       the cytosol to the nucleus when phosphorylated.
12874-CC       {ECO:0000250|UniProtKB:P48349}.

```

```
grep -A 6 -n -m 200 'CC   -!- SUBCELLULAR LOCATION.*(.*)'  /mnt/.../uniprot_sprot.dat

```

```
1560049:CC   -!- SUBCELLULAR LOCATION: Secreted. Note=Dense core vesicles (DCV).

2452012:CC   -!- SUBCELLULAR LOCATION: DNA-(apurinic or apyrimidinic site) lyase,
2452013-CC       mitochondrial: Mitochondrion. Note=Translocation from the
2452014-CC       cytoplasm to the mitochondria is mediated by ROS signaling and
2452015-CC       cleavage mediated by granzyme A. Tom20-dependent translocated
2452016-CC       mitochondrial APEX1 level is significantly increased after
2452017-CC       genotoxic stress (By similarity). The cleaved APEX2 is only
2452018-CC       detected in mitochondria. {ECO:0000250}.

9243083:CC   -!- SUBCELLULAR LOCATION: Isoform Soluble copper chaperone CopA(Z):
9243084-CC       Cytoplasm {ECO:0000250|UniProtKB:Q59385}.

38197671:CC   -!- SUBCELLULAR LOCATION: E(rns) glycoprotein: Host cell membrane
38197672-CC       {ECO:0000269|PubMed:28290554}; Peripheral membrane protein. Virion
38197673-CC       membrane {ECO:0000269|PubMed:28290554}; Peripheral membrane
38197674-CC       protein {ECO:0000305}. Note=The C-terminus membrane anchor of Erns
38197675-CC       represents an amphipathic helix embedded in plane into the
38197676-CC       membrane.

4924705:CC   -!- SUBCELLULAR LOCATION: Note=(Microbial infection) In the presence
4924706-CC       of M.tuberculosis EsxA-EsxB complex decreased amounts of B2M are
4924707-CC       found on the cell surface (PubMed:25356553).
4924708-CC       {ECO:0000269|PubMed:25356553}.

```

**02/18/2019**

For the DE line it is noted

> A block of DE lines may further contain multiple Includes: and/or Contains: sections and a separate field `Flags`: to indicate whether the protein sequence is a precursor or a fragment:

At present this seems to be related to the answer. When I have the file correctly parsed, then I can check for a correspondence as verification.

# Answer
> 1 votes
Thank you for asking this question. "Flag" is actually an obsolete notation which predates the introduction of proper evidence attribution in UniProtKB, and which escaped us in our efforts to update documentation to use the correct terminology. The correct term would now be "Evidence", and the next version of the user manual will have this corrected (release 2019\_06 due on July 3rd, 2019).

Before the current way of evidencing annotations was introduced in UniProtKB, the following "non-experimental qualifiers" or "flags" were used: By similarity, Potential, Probable. Absence of such a flag would generally imply that the information is based on an experiment.

When in doubt, please don't hesitate to contact the UniProt helpdesk.

---
Tags: file-formats, uniprot
---


# Question
Title: In Uniprot how do I know which variants are mutations and which are polymorphic?
In Uniprot how do I know which variants are mutations and which are polymorphic? I am looking at P38398 (BRCA1\_HUMAN) specifically.

Many thanks.

# Answer
> 2 votes
This answer is very late but may still help others: Have a look at the UniProt index of human polymorphisms and disease mutations at https://www.uniprot.org/docs/humsavar

Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.

---
Tags: uniprot
---


title: Fast processing of fastq data

text:
# Question
Title: Fast processing of fastq data
I am trying to write python script for customized filtering for fastq file (size \>3 GB). My proposed script is as follows:

```
def filtering(read):
    time.sleep(0.1)

    if len(read) >= 15 and \
       len(read) <= 30 and \
       np.mean(read.qual) >= 30 and \
       sum(i < 20 for i in read.qual) <=2:
           flag1 = True
           return  flag1

pool = multiprocessing.Pool(5)
pool.ncpus = 4
for file in os.listdir('/moRNA_data/trimmed_fastq'):
    if file.endswith('.fq'):
        fastq_file = HTSeq.FastqReader( os.path.join('/moRNA_data/trimmed_fastq',file))
        name = str(file.split('.')[0]) + '.fastq'
        file_out = open(os.path.join('/moRNA_data/filtered_fastq',name), "w" )
        for read in tqdm(fastq_file):
            flag1 = pool.apply_async(filtering,read)
            if flag1:
                read.write_to_fastq_file( file_out )
pool.close()
pool.join()

```

The problem with this script is that it works well with small file like upto 100-200 MB but when it goes to directory having multiple sample each having size more than 3 GB, it works well for first file but after that it hangs so severely that I had to restart the server \[ I can't even kill the process\]. I need suggestions to improve this script to make fastq file processing faster. Platform details is as follows:

```
OS : Ubuntu-18.02
Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              12
On-line CPU(s) list: 0-11
Thread(s) per core:  2
Core(s) per socket:  6
Socket(s):           1
NUMA node(s):        1
Vendor ID:           GenuineIntel
CPU family:          6
Model:               158
Model name:          Intel(R) Core(TM) i7-8700
CPU @ 3.20GHz
Stepping:            10
CPU MHz:             800.560
CPU max MHz:         4600.0000
CPU min MHz:         800.0000
BogoMIPS:            6384.00
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            12288K
NUMA node0 CPU(s):   0-11

```

Thanks.

# Answer
Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.

C++, Rust, Go, Nim and Julia can be as fast.

```
// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
// Compile with: gcc -O2 -o myprog this-file.c -lz
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
// STEP 1: declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r"); // STEP 2: open the file handler
    seq = kseq_init(fp); // STEP 3: initialize seq
    while (kseq_read(seq) >= 0) { // STEP 4: read sequence
        int i, sum, low;
        if (seq->qual.l == 0) continue; // no quality
        if (seq->seq.l < 15 || seq->seq.l > 30) continue;
        for (i = sum = low = 0; i < seq->seq.l; ++i) {
            int qual = seq->qual.s[i] - 33; // assuming Sanger quality
            sum += qual;
            if (qual < 20) ++low;
        }
        if (sum < seq->seq.l * 30.0 || low > 2) continue;
        printf("@%s\n", seq->name.s); // output can be made faster
        puts(seq->seq.s);
        printf("+\n");
        puts(seq->qual.s);
    }
    kseq_destroy(seq); // STEP 5: destroy seq
    gzclose(fp); // STEP 6: close the file handler
    return 0;
}

```

To apply the program to all `*.fq` in a directory (again, *not* tested):

```
mkdir -p filtered
ls *.fq | xargs -i echo ./myprog {} \> filtered/{} | sh

```

You can also replace the final `sh` with `parallel`, but the above should be fast enough.
> 6 votes

---
Tags: fastq, filtering, text-processing, data-preprocessing
---