# Question
Title: Convert Reactome Protein IDs to Pathway IDs?
I have a list of about 50k 'Protein IDs' from Reactome. Is there a simple way to get all the corresponding 'Pathway IDs' for each protein? What is the best service to use? (I'm guessing I can use the Reactome API, but I don't necessarily want to hit that 50k times...).

Protein IDs from reactome look like this:

```
R-HSA-49155
R-HSA-199420

```

The corresponding Reactome 'Pathway IDs' for those Protein IDs would be:

```
R-HSA-49155  R-HSA-110331
R-HSA-49155  R-HSA-110330
R-HSA-49155  R-HSA-110357
R-HSA-199420 R-HSA-1660499
R-HSA-199420 R-HSA-202424
R-HSA-199420 R-HSA-199418

```

# Answer
If you don't mind hitting it 50k times and are OK with python3...

```
from urllib import request
import json

def getPathways(proteinID):
    baseURL = 'http://reactome.org/ContentService/data/query'
    PathwayIDs = set()
    try:
        response = request.urlopen('{}/{}'.format(baseURL, proteinID)).read().decode()
        data = json.loads(response)
        if 'consumedByEvent' in data:
            for event in data['consumedByEvent']:
                PathwayIDs.add(event['stId'])
        if 'producedByEvent' in data:
            for event in data['producedByEvent']:
                PathwayIDs.add(event['stId'])
    except:
        pass
    return PathwayIDs

```

Usage would then be something like:

```
l = ['R-HSA-49155', 'R-HSA-199420', '']
for rid in l:
    ids = getPathways(rid)
    for _ in ids:
        print("{}\t{}".format(rid, _))

```

Which would produce:

```
R-HSA-49155 R-HSA-110239
R-HSA-49155 R-HSA-110240
R-HSA-49155 R-HSA-110238
R-HSA-49155 R-HSA-110356
R-HSA-199420    R-HSA-8948800
R-HSA-199420    R-HSA-6807106
R-HSA-199420    R-HSA-6807206
R-HSA-199420    R-HSA-6807126
R-HSA-199420    R-HSA-8847968
R-HSA-199420    R-HSA-8850997
R-HSA-199420    R-HSA-2321904
R-HSA-199420    R-HSA-8948775
R-HSA-199420    R-HSA-8944497
R-HSA-199420    R-HSA-6807134
R-HSA-199420    R-HSA-8850945
R-HSA-199420    R-HSA-8873946

```

Note that this will silently ignore invalid or missing IDs such as `''` (that's the `try` and `except` above. Note also that these are different pathway IDs than what you provided in your example. The main reason is that the protein IDs you showed are not always involved in the pathways IDs you showed (in my example, they always are).
> 6 votes

# Answer
On Reactome website they have at the download page mapping files (uniprot, ensembl, etc.), but unfortunately not for the protein IDs you are using (stable identifiers).

I had contact with their helpdesk, and they sent me a file containing all protein IDs to the pathways. Exactly what you need. I have asked them if they wanted to put it on their download page as well, but not sure if they want to do this. Meanwhile you can ask for the file as well, or get it from my google drive.

I assume you know how to get your Protein IDs and pathways from this file, using e.g., R?
> 3 votes

---
Tags: conversion, identifiers
---


# Question
Title: Is there a way to retrieve ENSEMBL IDs from a search query?
I would like to be able to batch download FASTA files from ENSEMBL. I normally would use the API to download them from the ENSEMBL gene IDs. However, the IDs aren't ones I know, only ones I would get from searching on the ENSEMBL website and collecting the IDs manually.

Is there a way I can use the API to query search terms and collect all results that are returned?

So for example, if I use the ENSEMBL website to search "HLA" I could collect the list of IDs (ENSG00000204252 etc. etc.) manually. HLA works okay using the API to query gene IDs because there's only a few terms (A, B, C plus non-classical for class I) but I was wondering if there was a way to directly access the search query one the homepage programmatically, as I may need to do it for some messier examples.

# Answer
> 1 votes
I would download the complete geneset FASTA file from ensembl and then parse out the genes you want with a script.

Alternatively, you could do similarly with API. Get all the gene names, select the ones you want, and then pull down each sequence.

---
Tags: database, reference-genome, ensembl, data-mining
---


# Question
Title: Building network for biomarker interaction analysis
I am looking for advice on how to best approach a problem I am faced with.

I have a dataset of numerous degradative biomarkers, clinical information and various other measures (from clinical trials). I would like to see how different biomarkers relate to each other, and effect their levels (high/low etc).

My goal is to gain an understanding of how these biomarker levels fluctuate in different disease profiles - and understand if there is a relationship between each of them (co/independent). The limitation is that I do not have access to genetic or environmental data so this will be a "random effect". There are around 1000 marker samples and 10's of biomarkers mixed between clinical and biochemical. There are also numerous time points to be incorporated.

I have thought about implementing something similar to a protein-protein interaction network, although this will of course not involve physical interactions, but more pathway interactions for further hypothesis driven research.

I wanted to firstly: validate if this project actually makes sense and secondly: ask for some advice around how to best implement this.

I have experience implementing neural/deep nets, directed and undirected networks (Bayes) in python, but not sure exactly how to go about this.

Thanks!

# Answer
> 2 votes
You can start by looking at the correlation (if data is non-parametric) of each variable with each variable of each data point you have. One thing that might change is the correlation between variables along the disease.

Next you can determine which relationships are constant along the disease progression/clinical trial. Identify also the ones that don't show a constant relationship. To do so you can use WGCNA (although with just 10 variables is too low to expect a scale free network you can use to correlate the variables and group them by correlation profile)

Then the question is which variables are more important, those who kept the same relationship or those which change? Use the eigengenes of those variables (separately by those that have a constant relationship and those which don't) to estimate the dependency with the parametric variables. Or you simply can do an ANOVA of all the variables with each non-parametric variable.

---
Tags: proteins, networks, pathway, interactions
---


# Question
Title: How can I find mutations associated with disease in human histone residues?
I would like to look if there are mutations in residues of human histones associated with any disease. For instance, if a mutation in residue K6 (lysine 6) of histone H2A1A is associated with any human disease, (with PUBMED evidence, preferably), but in a systematic way (for all types of histones and all residues).

The purpose of this is to get some candidate residues to study. I have checked the OMIM database but I didn't find anything relevant. It seems that there is no systematic way to connect specific residues to diseases and search for them.

Is there any web database or API available for this purpose?

# Answer
> 5 votes
**Ensembl** contains this information: When you go to the “phenotype” menu item of a given gene, you will see a list of variants (potentially after clicking on “ALL associated variants”) with their associated genomic position.

You can subtract the transcript’s start position from that position to find out the residue (of the translated sequence, and be careful to subtract the lengths of all preceding introns!). Unfortunately this seems to require some manual work.

Alternatively you can query the individual SNPs, which will tell you their position in all transcripts (in this case: Leu 298/68/244):

# Answer
> 3 votes
I found DisGeNET useful for my purpose, a database that associates genes with diseases and, if known, gene variants with diseases. It integrates several sources of information for computing a score of association of gene/variation with the disease ( Piñero et al., 2016, Piñero et al., 2015):

* **Gene Disease Associations** from UniProt, PsyGeNET, ClinVar, Orphanet, the GWAS Catalog, CTD (human data), and Human Phenotype Ontology, RGD, MGD, and CTD (mouse and rat data) GAD, LHGDN and BeFree
* **Variants Disease Associations** from ClinVar, the GWAS Catalog, Uniprot, GAD, and from BEFREE

If available, it provides the Pubmed ID of articles supporting the proposed association and the rs id from the dbSNP for variants associated to disease.

In addition it provides an experimental R package, disgenet2r, that allows to query the database in a programmatic manner and that provides some plot functions to generate representations of the associations founded.

Using this approach, I have been able to found several human histone-disease associations (not my original purpose), together with some concrete histone residues that are mutated in specific diseases.

---
Tags: database, variation, api
---


# Question
Title: Profile of conservation between bacterial sequences and human protein
My goal is to make a conservation plot between bacterial sequences and a human protein.

So far, I have a FASTA file of the protein, and a FASTA file with the sequences of the proteins from the BLAST results.

In my initial attempts, I have tried to do this using msa software:

I have downloaded clustalx and run a profile alignment using the single protein sequence as Profile 1 and the sequences from the BLAST results as Profile 2, and selected the option "Align Sequences to Profile 1".

I'm having difficulty interpreting and dealing with the results. Since the majority of the sequences aren't aligned, the results are a total mess, with a huge number of dashes added.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

ie:

I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.

Is there a better way to achieve this or to salvage the results from the profile alignment?

Thanks.

# Answer
> 2 votes
I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the *same genome* so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.

You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

> My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.

If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py

Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.

For example, given this MSA:

```
    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
           ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
           ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
           EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
           DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH

```

You'd get this plot:

> I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.

What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.

As you say in the comments,

> My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

1. Blast sequence of interest.
2. Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
3. Align all the sequences.
4. Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.

---
Tags: sequence-alignment
---


# Question
Title: how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation?
For DNA/RNA quantification machines like the Bioanalyzer or TapeStation, the DNA Integrity Number (DIN) or RNA Integrity Number (RIN) numbers are quoted as a measure of the fragmentation of the material.

How is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation? Is there any open-source software to parse/collate the data from the machines?

# Answer
> 2 votes
Note that DIN (DNA Integrity Number) and RIN (RNA Integrity Number) are different scores. I've had trouble in the past finding code or formulas for calculating both of these scores, which is very frustrating considering how frequently they are used in research papers and for NGS QC. The closest that I have been able to get to this are demonstrative graphs for various RIN numbers and DIN numbers, with a description of the variables that are included in the model.

The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.

Just in case it's helpful, here is an extract of that code, only including lines that seem to be directly involved in calculating the actual B\_RIN values. This is not runnable, because it depends on particular input files, and it's still pretty unintelligible:

```
myfiles<-list.files()
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
    ##mod chip numb
    targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
    str3<-paste0(substr(targets2[1,2],1,67),"_S")
    samp_num_chip<-length(grep(str3,myfiles))
    ##mod sample numb
    str1<-substr(targets2[1,2],1,67)
    for(j in 1:samp_num_chip)
    {
        iii<-j
        str5<-paste0(str1,"_Sample",iii,".csv")
        samples1[j,i]<-str5
    }
}
samples1<-samples1[grep("2100",samples1)]
##RIN Calculator#######refined peakfinder####################
##read bioanalyzer data into a matrix called dta
##since the total RNA and mRNA assays run differently
##skip more lines for the total RNA assay
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
    x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
    dta[,i] <- x[,2]
}
time<- x[,1]
for(samples in 1:length(samples1))
{
    ##store max fluor in max.peaks matrix
    lad<-samples
    ti1<-17
    ti2<-65
    p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
    peaks <- matrix(NA, nrow=p.row, ncol=1)
    max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
    for (i in 1:length(seq(ti1,ti2,.05)))
    {
        ii<-seq(ti1,ti2,.05)[i]
        x0<- ii
        x0<-round(x0,digits=3)
        y0<- ii
        x<- ii+.05
        x<-round(x,digits=3)
        y<- ii
        max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
    }
    ##human or mouse
    cent<-max(max.peaks)*.6
    ##flatworm
    ##cent<-max(max.peaks)*.05
    ##define the middle of the segments
    max.avg<-seq(ti1+.025,ti2+.025,.05)
    ##find time associated with peaks > 5 fu and fill seconds matrix with them
    ##12.5% of max fu incase chip is fu
    sizes <- which(max.peaks>cent)
    gt35<-which(max.avg[sizes]>35)
    lt55<-which(max.avg[sizes]<55)
    tfr<-which(duplicated(c(gt35,lt55))==TRUE)
    area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],
                max.peaks[sizes][c(gt35,lt55)[tfr]],
                max.peaks[sizes+1][c(gt35,lt55)[tfr]])
    time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],
                 max.avg[sizes][c(gt35,lt55)[tfr]],
                 max.avg[sizes+1][c(gt35,lt55)[tfr]])
    t1<-time2[which(area[,1]<cent),1][1]
    t2<-time2[which(area[,1]<cent),1][2]
    t3<-time2[which(area[,3]<cent),3][1]
    t4<-time2[which(area[,3]<cent),3][2]
    s18a<-ifelse(is.na(t1,-30,round(t1,digits=1)))
    s18b<-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
    s28a<-ifelse(is.na(t2,-30,round(t2,digits=1)))
    s28b<-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
    eta<-which(time==s18a)
    etb<-which(time==s18b)
    tea<-which(time==s28a)
    teb<-which(time==s28b)
    s18<-sum(dta[eta:etb,lad])
    s28<-sum(dta[tea:teb,lad])
    qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
}

```

Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.

# Answer
> 1 votes
The best I could find by looking at github.com is:

https://github.com/brianfleharty/Agilent-BioAnalyzer-RIN-R-code-for-Fly/blob/master/RNA\_RIN.R

---

```
#set the working directory
dir<- "H:\\B_RIN\\emd"
setwd(dir)
myfiles<-list.files()
which(substr(myfiles,69,79)=="Results.csv")
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
length(n)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str1
str2<-"_S"
str3<-paste(str1,str2,sep="")
str3
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
#mod chip numb
targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str2<-"_S"
str3<-paste(str1,str2,sep="")
samp_num_chip<-length(grep(str3,myfiles))
#mod sample numb
str1<-substr(targets2[1,2],1,67)
str4<-"_Sample"
for(j in 1:samp_num_chip)
{
iii<-j
end.st<-".csv"
str5<-paste(str1,str4,iii,end.st,sep="")
samples1[j,i]<-str5
}
}
samples1
grep("2100",samples1)
samples1<-samples1[grep("2100",samples1)]
samples1
write.table(samples1,file=paste(dir,"\\ba_lane.txt",sep=""), sep="\t",col.names=NA)
#RIN Calculator#######refined peakfinder####################
#read bioanalyzer data into a matrix called dta
#since the total RNA and mRNA assays run differently
#skip more lines for the total RNA assay
ramp <- colorRamp(c("darkmagenta","white"))
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
dta[,i] <- x[,2]
}
time<- x[,1]
#plot ladder and define peak threshold in HD
par(mfrow=c(3,4))
for(samples in 1:length(samples1))
#for(samples in 1:nrow(elec))
{
#which(dta[,lad] > 5)
#draw some segments and store max fluor in max.peaks matrix
lad<-samples
ti1<-17
ti2<-65
p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
peaks <- matrix(NA, nrow=p.row, ncol=1)
max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
for (i in 1:length(seq(ti1,ti2,.05)))
{
ii<-seq(ti1,ti2,.05)[i]
x0<- ii
x0<-round(x0,digits=3)
y0<- ii
x<- ii+.05
x<-round(x,digits=3)
y<- ii
#segments(x0,y0,x,y)
max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
}
#human or mouse
cent<-max(max.peaks)*.6
#flatworm
#cent<-max(max.peaks)*.05
plot(time,dta[,lad],type="h",col="dodgerblue",xlab="seconds",ylab="Fluorescence Units",xlim=c(17,75),ylim=c(0,(max(max.peaks)*1.2)),main=elec[lad,1])
abline(h=cent,col="grey")
#define the middle of the segments
max.avg<-seq(ti1+.025,ti2+.025,.05)
#plot the maximum value of each segment
max.peaks
lines(max.avg,max.peaks,type="b",col="deepskyblue")
#lines(max.avg,min.peaks,type="b",col="blue")
#find time associated with peaks > 5 fu and fill seconds matrix with them
#12.5% of max fu incase chip is fu
p1<-which(max.avg==38)
p2<-which(max.avg==55)
sizes <- which(max.peaks>cent)
max.avg[sizes]
gt35<-which(max.avg[sizes]>35)
lt55<-which(max.avg[sizes]<55)
#as.numeric(duplicated(c(gt35,lt55)))
tfr<-which(duplicated(c(gt35,lt55))==TRUE)
c(gt35,lt55)[tfr]
max.avg[sizes][c(gt35,lt55)[tfr]]
max.peaks[sizes][c(gt35,lt55)[tfr]]
points(max.avg[sizes+1][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes-1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],col="deeppink3",pch=19)
#colnms2<-c("max.avg[sizes]","max.peaks[sizes]","max.avg[sizes+1]","max.peaks[sizes+1]","max.avg[sizes-1]")
max.avg[sizes]
max.peaks[sizes]
max.avg[sizes+1]
max.peaks[sizes+1]
max.avg[sizes-1]
max.peaks[sizes-1]
area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]])
time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.avg[sizes][c(gt35,lt55)[tfr]],max.avg[sizes+1][c(gt35,lt55)[tfr]])
p1<-area[which(area[,1]<cent),1][1]
t1<-time2[which(area[,1]<cent),1][1]
p2<-area[which(area[,1]<cent),1][2]
t2<-time2[which(area[,1]<cent),1][2]
#area[which(area[,2]<cent),2]
#area[which(area[,2]<cent),2]
p3<-area[which(area[,3]<cent),3][1]
t3<-time2[which(area[,3]<cent),3][1]
p4<-area[which(area[,3]<cent),3][2]
t4<-time2[which(area[,3]<cent),3][2]
#abline(v=area[which(area[,3]<cent),3])
#abline(h=11)
#abline(v=42.75)
v1<-t1
v2<-t2
v3<-t3
v4<-t4
abline(v=v1)
abline(v=v2)
abline(v=v3)
abline(v=v4)
#time[area[which(area[,1]<cent),1][1]==dta[,lad]]
#time[area[which(area[,3]<cent),3][1]==dta[,lad]]
#time[area[which(area[,1]<cent),1][2]==dta[,lad]]
#time[area[which(area[,3]<cent),3][2]==dta[,lad]]
#abline(v=time[area[which(area[,1]<cent),1][1]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][1]==dta[,lad]])
#abline(v=time[area[which(area[,1]<cent),1][2]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][2]==dta[,lad]])
#abline(h=area[which(area[,3]<cent),3][1])
#abline(h=area[which(area[,3]<cent),3][2]+3)
s18a<-round(t1,digits=1)
s18b<-round(t3,digits=1)
s28a<-round(t2,digits=1)
s28b<-round(t4,digits=1)
if(is.na(s18a)){s18a<-30}
if(is.na(s18b)){s18b<-35.05}
if(is.na(s28a)){s28a<-30}
if(is.na(s28b)){s28b<-30.05}
eta<-which(time==s18a)
etb<-which(time==s18b)
tea<-which(time==s28a)
teb<-which(time==s28b)
#dta[eta:etb,lad]
#dta[tea:teb,lad]
s18<-sum(dta[eta:etb,lad])
s28<-sum(dta[tea:teb,lad])
s28/s18
qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
text(60,cent*.75,c(format(qc.mat[samples],digits=2),"\n\nB-RIN"))
}
qc.mat
###end of ratio calc
###
###
#barplot(qc.mat,beside=T,ylim=c(0,11))
#set the working directory
#results <- read.delim("results-files.txt", sep=",",header=F,skip=0,as.is=T)
#targets <- read.delim(results[1,], sep=",",header=F,skip=14,as.is=T)
#targets[1:5,]
tbl<-which(targets[,1]=="Overall Results:")
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
#read in the targets file
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
mat1 <- matrix(NA, ncol=length(colnms), nrow=length(samples)*length(results))
colnames(mat1)<-colnms
rownames(mat1)<-rep(samples,length(results))
for(i in 1:length(results))
{
targets <- read.delim(results[i], sep=",",header=F,skip=14,as.is=T)
targets[1:5,]
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
tbl<-which(targets[,1]=="Overall Results:")
length(which(targets[,1]=="Overall Results:"))
tbl[1]
tbl[1]
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
for (j in 1:length(samples))
{
mat1[((length(samples)*(i-1))+j),]<-as.matrix(targets[(tbl[j]+1):(tbl[j]+4),2])
}
}
mat1
cbind(mat1,qc.mat)
colnames(qc.mat)<-"B-RIN"
write.table(cbind(mat1,qc.mat),file=paste(dir,"\\RIN_VS_B-RIN.txt",sep=""), sep="\t",col.names=NA)

```

---
Tags: phylogenetics, quantification
---


# Question
Title: What tools can I use for a bacterial core/pan genome pipeline?
I want to perform a genome comparison on a group of isolates. I want to look into two broad groups of taxa and compare the accessory genome in each group. I have been using prokka (v1.12) and roary (v3.8.2) to do this but it appears the accessory\_binary\_genes.fa file is actually an untrue representation.

**Note:** gene\_presence\_absence.Rtab does contain all the full presence/ absence for accessory gene sets. Despite this Im still unhappy with the nomenclature of the gene groups \[issue for another day\]

> \[github issue 335\] Your best to ignore the accessory\_binary\_genes.fa file. It is just for creating a quick and dirty tree with FastTree. The file itself is filtered to remove very common and not common variation to speedup the tree generation, hence the difference in numbers.
>
> The top 5% and bottom 5% are excluded. It is truncated at 4000 genes.

I've been looking into alternative pipelines and a new software BPGA looks promising. Does anyone have experience with this?

I essentially want a tool which will give me the core and accessory gene sets, without the noise from partial gene hits.

# Answer
> 3 votes
LS-BSR should be able to give you what you're looking for. See the article.

After you run the primary analysis, there is a simple documented workflow for splitting the pangenome into core and accessory, based on a user defined threshold. I'm the developer, so can help if you run into any problems.

# Answer
> 2 votes
Roary also takes into account paralogs, so sometimes two core genes are split into different groups based on their neighbour genes and they end up with different nomenclature (group\_\*...). As suggested by Andrew Page in the github issue I would consider the gene\_presence\_absence.Rtab (this contains all the orthologous genes) and remove rows corresponding to vectors only containing 1s (core genes). In this way you will have a matrix of 1 and 0 corresponding to presence/absence of a particular accessory gene in your isolates.

---
Tags: genome, software-recommendation
---


# Question
Title: What methods should I use from PythonCyc API to query metabolites in BioCyc database?
I am using **PythonCyc** API in order to write a query for metabolites in **BioCyc**. The purpose of this API is to communicate with the database software of BioCyc- Pathway Tools. Pathway Tools is in lisp therefore, PythonCyc creates a bridge between python and common lisp. To use this this API one must first create a PGDB object with a specified organism ID (orgid). In the example below I create a PGDB with orgid "meta". After this, I am able to call methods from PythonCyc with this object:

```
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('sucrose')
print pathways

```

The query for the metabolite above, 'sucrose', provides a list a pathways:

```
[u'|PWY-7347|', u'|SUCSYN-PWY|', u'|PWY-7238|', u'|PWY-6527|', u'|PWY-5343|', u'|PWY-3801|', u'|PWY-7345|', u'|PWY-862|', u'|SUCUTIL-PWY|', u'|PWY66-373|', u'|PWY-621|', u'|PWY-822|', u'|SUCROSEUTIL2-PWY|', u'|PWY-6525|', u'|PWY-6524|', u'|PWY-5337|', u'|PWY-5384|']

```

However if I switch this metabolite name to a common amino acid name, such as 'valine':

```
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('valine')
print pathways

```

The query provides an error message that states "non coercible frame", which means it cannot find this ID. This is a common metabolite, just as sucrose, which should very well have pathway entries in database, however there are none found with this method. I have also tried synonyms of valine and other methods such as "all\_pathways(cpds)" stated in API, which give me the same error message. What is the reason for this? Does BioCyc not list any pathways for amino acids? Or am I using an incorrect method to access amino acid information?

Link to PythonCyc API

# Answer
> 3 votes
**UPD**: As suggested by Llopis, maybe the correct identifier is VAL?

I think there is no metabolite with such name in the database. I tried their search, the correct name seems to be `L-valine`. I expect that they use similar naming for other amino acids. If you cannot find something, you can use their compound search through the website.

# Answer
> 3 votes
Pathway-tools has a unique identifier (called a Frame-ID) for each metabolite in its database. You can get the list of all compounds with this api call:

`all_cpds = meta.get_class_all_instances('|Compounds|')`

for each compound, you can also get its common name, or all its synonyms:

`cpd_names = {}
for cpd_frame_id in all_cpds:
 cpd_names[cpd_frame_id] = meta.get_name_string(cpd_frame_id)`

---
Tags: python, api
---


# Question
Title: Obtaining all protein sequences with a particular domain architecture from Pfam
I want to get the alignment of chain A of 1kf6 (PDB ID) from the pfam database here. This protein chain has two main domains (FAD\_binding\_2 and Succ\_DH\_flav\_C). In pfam there is a link to one of these domains and after clicking in one of the mentioned domains in the table below the page, another page comes at the top of which there is a link to architectures.

For example, if we click on "FAD\_binding\_2" in the table, there are 220 architectures that have this domain and after clicking on that button (at the top of page), a page opens here. In this page 7471 sequences with the following architecture: FAD\_binding\_2, Succ\_DH\_flav\_C, resembles most to my chain (because it has both FAD\_binding\_2 and Succ\_DH\_flav\_C domains).

If we want to see all 7471 sequences with that architecture we can click on the show button under that architecture and finally we can see all of that 7471 architecture resembling our chain A sequence. Down the page we can see other architectures that all have FAD\_binding\_2, but some of them do not have both FAD\_binding\_2 and Succ\_DH\_flav\_C domains (or some of them have other extra domains that we are not interested in) and so we do not need their sequences and their alignments.

When we want to get the alignment of full sequences at the left side of the page, we can see that we can just get the alignment of all 15696 sequences and it seems that there is no way to get the alignment of our favourite 7471 sequences. I would like to know if there is any way to get the alignment or at least the fasta file of our favourite 7471 sequences (which have both of our favourite domains not just one of them)?

# Answer
I believe you can do this using SQL files of the database.

First you need to get the integer accession of this architecture. `architecture` table has the following columns:

* `auto_architecture` (integer id, **this is what we need**)
* `architecture` (a string describing the architecture)
* `type_example` (example sequence)
* `no_seqs` (number of sequences)
* `architecture_acc` (accessions of the architectures)

You need the following line:

```
3900719357      FAD_binding_2~Succ_DH_flav_C    Z9JRB3  7471    PF00890 PF02910

```

Now you need to get the sequences with this architecture id. First, download the file with all the sequences (**warning**, this one is more than 7Gb). With this one you need column #1 (pfam sequence accesion), #12 (sequence itself) and #16 (architecture integer id).

If you do not want to import files in the database, you can use `awk` to export all the sequences to `fasta` format.

```
zcat pfamseq.txt.gz | awk -F\\t '{ if ($16=="3900719357") print ">" $1 "\n" $12}' > sequences.fasta

```

Of couse, if you need to perform multiple queries like this, it is much easier to import tables into the database. Then you can perform both queries in a single SQL operation (assuming you know the architecture string):

```
SELECT pfamseq_acc, sequence
    FROM architecture
    JOIN pfamseq ON pfamseq.auto_architecture=architecture.auto_architecture
    WHERE architecture="FAD_binding_2~Succ_DH_flav_C";

```

(I haven't tried this SQL query, so maybe you need to correct it a bit)
> 2 votes

# Answer
An alternative way is using the Hmmer website, and do a hmmsearch with the Pfam accession, forex. (PF02910). Navigate to the domain tab and find the domain architecture that you are interested in and press view scores. That way, only the sequences you are interested in will be downloadable under the tab downloads.

Step-by-Step guide:

1- Go to http://www.ebi.ac.uk/Tools/hmmer/

2- Search with hmmsearch http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch

3- Use Accession search and type in PF02910

4- Navigate to Domain

5- To the right side of 9742 sequences with domain architecture: FAD\_binding\_2, Succ\_DH\_flav\_C

Press View scores.

6- It should show (Your results have been filtered) and Query matches (9742)

7- Navigate to Download and download the sequences in the format of interest.

8- Repeat for domain architectures of interest.

I hope you find this helpful.

Kind Regards,

Sara
> 3 votes

# Answer
That does not appear to be possible using PFAM. I would suggest trying a different database. I managed to perform the same query using SMART (http://smart.embl-heidelberg.de/), via entering "Pfam:FAD\_binding\_2 AND Pfam:Succ\_DH\_flav\_C" under "Domain selection" in the "Architecture Analysis" section of the search screen, which resulted in 9921 protein sequences with those two domains as you can see here.

On the results screen you can select all the sequences and under "Action" at the top of the screen choose "Download Protein Sequences as FASTA file".
> 2 votes

---
Tags: sequence-homology, domains
---


# Question
Title: Calculating the charge of a peptide computationally
I was wondering how I can calculate the charge of a protein peptide (e.g. "RKTTLVPNTQTASPR") computationally in R or another tool.

# Answer
A quick google search turns up protcalc, which is able to give a nice pH-dependent table of peptide charges (yours ranges from 3.1 at pH 4 to 1.5 at pH 10). It's on sourceforge, so hopefully the source code (and maybe a publication) is somewhere in there. Granted, this isn't R, but it's not clear from your post how necessary that really is.

BTW, in general "best and most accurate" turns out to be highly subjective (at least the "best" part of that).
> 3 votes

# Answer
Another quick Google search reveals the Peptides R package.

```
library(Peptides)
s <- "RKTTLVPNTQTASPR"
charge(s)

[1] 2.997683

```

Optional arguments to `charge()` are `pH` (default = 7) and `pKscale` (choice of nine, default = "Lehninger"). See `?charge` for details.

Another (non-R) option is the EMBOSS suite program iep. I highly recommend installing EMBOSS, it provides a suite of useful command-line tools with output that can be easily read and processed by other sofware, including R.
> 3 votes

# Answer
I've used PropKa in the past to get isoelectric points. Pretty simple to use:

https://github.com/jensengroup/propka-3.1
> 1 votes

---
Tags: r, proteins
---


# Question
Title: Does Prokka or gff3 change the sequence / subject ID of fasta files?
I used prokka to create a database for some fasta files but I noticed a strange difference between the prokka and fasta files. Normally the seqid for the gff3 output for prokka has seqid's that match the sequence id's of the fasta files. The files I just produced with prokka have a different seqid than my fasta files. My fasta files look like this

```
>NODE_108_length_645_cov_0.679537_ID_215

```

However the prokka output gff3 file for the same fasta file uses a seqid like this

```
gnl|X|JHIOGLGG_28

```

What would cause this difference? And do the last digits still correlate to the same seqid's even though the string before it has changed?

# Answer
> 3 votes
Prokka is an annotation tool. It takes your contigs file with fasta headers `>NODE_108...`and searchs for ORFs and their putative functions. Prokka produces submission ready output, the new fasta headers you see `>JHIOGLGG_28...`are locus tags and gene description from the annotation and are not related to contig node data. Please see the github information page for more info.

---
Tags: fasta, gff3, sequence-analysis
---


title: What does the "X_OS_IND" column mean?

text:
# Question
Title: What does the "X_OS_IND" column mean?
I want to survival analysis using the subset of TCGA LUAD dataset, which identifier are located here.

```
> head(LUAD_clinicalMatrix[,c("X_EVENT","X_OS_IND")])
  X_EVENT X_OS_IND
1      NA       NA
2       0        0
3       1        1
4       0        0
5       0        0
6       0        0

> all(LUAD_clinicalMatrix$X_EVENT==LUAD_clinicalMatrix$X_OS_UNIT)
[1] FALSE

```

I want to know which identifier means the sample is alive or dead, and "1" means alive or dead?

# Answer
According to this link, 1=death 0=censor and null=no data.
> 5 votes

---
Tags: r, file-formats
---


# Question
Title: Where can I find a GWAS data set about genotypes and phenotypes for Zebrafish and rice?
My goal is to work on GWAS method development as a research project. However, human genome data are usually confidential because of the identification problem, so it's very hard to get them. (application takes a long time and I may not get it at the end.)

Therefore, maybe I can try to play with the data of other species first. So I wonder where I can find some GWAS data for Zebrafish and rice. (The reason I am asking for these two is that I find researchers tend to perform genome analysis to mice, Arabidopsis Thaliana, zebrafish, and rice, and I already found where to download the GWAS data for the first two.)

I can find several data sets (like ZFIN) that have the annotated genome information, but that's only one genome. Where can I download the raw data (like genotypes of 500 zebrafishes and corresponding phenotypes) to start my own analysis?

# Answer
> 1 votes
I think it will be more productive for you if you choose one of the public human datasets and practise your GWAS methods on that. A couple of examples examples:

* 1000 genomes -- use "population" as your trait of association. You'll find that you get association across the whole genome, so the analysis method would change from "what bits are associated" to "what are the bits that are most consistently associated"
* NCI Genomic Data Commons -- open access (and controlled access) data associated with cancer

---
Tags: database, gwas
---


# Question
Title: ribo interval list for picard ensembl mm9 gtf
Could some one please help me with understanding how to generate the ribosomal interval list that is required to when using picard metrics this step

```
java -jar picard.jar CollectRnaSeqMetrics \
      I=input.bam \
      O=output.RNA_Metrics \
      REF_FLAT=ref_flat.txt \
      STRAND=SECOND_READ_TRANSCRIPTION_STRAND \
      RIBOSOMAL_INTERVALS=ribosomal.interval_list

```

I have a couple of bam files that were aligned using mm9 ensembl gtf as annotation, i did a lot of search and found this link Ribosomal Intervals For Collectrnaseqmetrics but this link or any other source does not explain how to produce or find ribo interval list from ensembl annotation

# Answer
At least with Ensembl annotations, you can `grep rRNA foo.gtf > rRNA.gtf` to get the annotated rRNAs. Note however, that the 45S cassette isn't in the mouse reference genome, so you'll never have alignments to 18S, 5.8S, or 28S rRNA. I always get the feeling that the Collect\*Metrics programs are made with human sequences in mind, since the cassette is included in the human reference genome.

Anyway, while you can use the annotation and reference sequence to make picard happy, make sure to append the Rn45S sequence to your genome if you want to actually measure rRNA levels. In my experience this is only useful for rRNA-depleted samples, since the rates can vary wildly there.

**Edit**: To make a file compatible with picard's "interval list", one can do the following with the GTF file:

```
grep rRNA foo.gtf | awk 'BEGIN{OFS="\t"}{print $1, $4-1, $5}' > foo.bed

```

According to the picard documentation a BED file can be used as input and the above `awk` code extracts the appropriate columns and switches to 0-based coordinates used in the BED format.
> 3 votes

---
Tags: rna-seq
---


# Question
Title: Counts obtained by featureCounts seem much less than observed coverage
I have surprisingly low counts when running `featureCounts` on some (single-end) RNA-seq data mapped on *C. elegans* genome using `hisat2`.

To more easily show the problem, I generated a small subset of the bam file and of the annotation file I'm using. Here is what I can see when loading these two files on IGV:

The coverage for the "his-11" and "his-15" genes peak at around 3600. For "his-12" and "his-16", it peaks at more than 1000

Here are the content of my annotation file:

```
$ cat his_11-16.gtf
II  ensembl transcript  13817759    13818143    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001890"; gene_name "his-16"; gene_source "ensembl"; gene_version "1"; p_id "P8354"; transcript_biotype "protein_coding"; transcript_id "ZK131.10"; transcript_name "ZK131.10"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS11539";
II  ensembl transcript  13818371    13818739    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001889"; gene_name "his-15"; gene_source "ensembl"; gene_version "1"; p_id "P8185"; transcript_biotype "protein_coding"; transcript_id "ZK131.9"; transcript_name "ZK131.9"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52436";
II  ensembl transcript  13819626    13819937    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001888"; gene_name "his-14"; gene_source "ensembl"; gene_version "1"; p_id "P21890"; transcript_biotype "protein_coding"; transcript_id "ZK131.8"; transcript_name "ZK131.8"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS50308";
II  ensembl transcript  13820210    13820620    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001887"; gene_name "his-13"; gene_source "ensembl"; gene_version "1"; p_id "P2149"; transcript_biotype "protein_coding"; transcript_id "ZK131.7"; transcript_name "ZK131.7"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS35106";
II  ensembl transcript  13821198    13821581    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001886"; gene_name "his-12"; gene_source "ensembl"; gene_version "1"; p_id "P26082"; transcript_biotype "protein_coding"; transcript_id "ZK131.6"; transcript_name "ZK131.6"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS23693";
II  ensembl transcript  13821778    13822314    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001885"; gene_name "his-11"; gene_source "ensembl"; gene_version "1"; p_id "P10535"; transcript_biotype "protein_coding"; transcript_id "ZK131.5"; transcript_name "ZK131.5"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52792";

```

The `featureCounts` run:

```
$ featureCounts -a his_11-16.gtf -o his_11-16_counts.txt -t transcript -g gene_name -O his_11-16_sorted.bam

        ==========     _____ _    _ ____  _____  ______          _____
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
      v1.5.2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S his_11-16_sorted.bam                           ||
||                                                                            ||
||             Output file : his_11-16_counts.txt                             ||
||                 Summary : his_11-16_counts.txt.summary                     ||
||              Annotation : his_11-16.gtf (GTF)                              ||
||      Dir for temp files : ./                                               ||
||                                                                            ||
||                 Threads : 1                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : counted                                          ||
||   Min overlapping bases : 1                                                ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file his_11-16.gtf ...                                     ||
||    Features : 6                                                            ||
||    Meta-features : 6                                                       ||
||    Chromosomes/contigs : 1                                                 ||
||                                                                            ||
|| Process BAM file his_11-16_sorted.bam...                                   ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 35037                                                     ||
||    Successfully assigned reads : 1849 (5.3%)                               ||
||    Running time : 0.00 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
|| Summary of counting results can be found in file "his_11-16_counts.txt"    ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

```

And the resulting counts file:

```
$ cat his_11-16_counts.txt
# Program:featureCounts v1.5.2; Command:"featureCounts" "-a" "his_11-16.gtf" "-o" "his_11-16_counts.txt" "-t" "transcript" "-g" "gene_name" "-O" "his_11-16_sorted.bam"
Geneid  Chr Start   End Strand  Length  his_11-16_sorted.bam
his-16  II  13817759    13818143    -   385 869
his-15  II  13818371    13818739    +   369 3
his-14  II  13819626    13819937    -   312 953
his-13  II  13820210    13820620    +   411 23
his-12  II  13821198    13821581    -   384 423
his-11  II  13821778    13822314    +   537 1

```

**Why such a low assignment rate, why so little counts?**

I see that `hisat2` decided to split an important proportion of the reads between orthologs, like "his-16" and "his-12", but my understanding is that the default values for the `minOverlap` and `fracOverlap` parameters should ensure that such split reads will be counted:

```
  --minOverlap <int>  Minimum number of overlapping bases in a read that is
                      required for read assignment. 1 by default. Number of
                      overlapping bases is counted from both reads if paired
                      end. If a negative value is provided, then a gap of up
                      to specified size will be allowed between read and the
                      feature that the read is assigned to.

  --fracOverlap <float> Minimum fraction of overlapping bases in a read that is
                      required for read assignment. Value should be within range
                      [0,1]. 0 by default. Number of overlapping bases is
                      counted from both reads if paired end. Both this option
                      and '--minOverlap' option need to be satisfied for read
                      assignment.

```

What did I get wrong?

# Answer
> 4 votes
Given the high level of multimapping in this region, you'll need to use the `-M --primary` options if you want to keep many of the alignments. I would be very hesitant to use these numbers as input for DESeq2 or similar programs, since it's fairly questionable whether one should fully trust the "randomness" of the aligner's assignments. I'm more comfortable using something like salmon for such cases.

---
Tags: read-mapping, coverage
---


title: insert size pre and post trimming

text:
# Question
Title: insert size pre and post trimming
I have a problem here with my rna seq data:

# Sequencing details:

rRNA was removed, followed by cDNA preparation and generation of stranded libraries using the TruSeq Stranded Total RNA Sample Prep Kit. Sequencing performed on the HiSeq2500 platform (Illumina) to generate 2 × 125 bp paired-end reads

# Alignment and preprocessing

reads were aligned using tophat2(std aligner in pipeline at core (`--library-type fr-firststrand`), unfortunately the original unaligned files were purged and i only have access to this aligned file which i converted to fastq using following steps

1. `samtools sort -n sample.bam -o sample_sorted.bam`
2. `bedtools bamtofastq -i sample_sorted.bam -fq sample_1.fq -fq2 sample_2.fq` (get a lot of mate skipping errors here)
3. check for these reads for adapters and trim them and again use hisat2 (with `--rna-strandness RF`)

# trimmed reads

Reads were trimmed by passing parameters `--adapters adapters.file --adapter-trim-end RIGHT --length-dist --threads 12 --adapter-min-overlap 7 --max-uncalled 250 --min-read-length 25 --` to FLEXBAR version 2.4

## adapters for all samples using multiqc

## fastqc sequence length pre trimming

## fastqc sequence length post trimming

## rmats error

> Incorrect readLength. sample.bam has a read length of 114, while readLength param is 125

## picard error on aligned bam from hisat2)

> ProcessExecutor  "Not creating insert size PDF as there are duplicated header names: All\_Reads"

## snapshot of picard insertsize metrics file

# Questions

1. Is my approach flawed at some step?
2. why do i get different read lengths error for each sample when i process it through rmats when i specify len as 125?
3. why is no histogram being produced when using picard insertsize metrics?
4. is normal to see the graph shift to the right after trimming?

# insertsize for all samples (data from picard collectinsertsize compiled using multiqc)

# Answer
1. It's not an unreasonable approach. rMATs is rather picky about its input, but it seems you've noticed that.
2. rMATs can't handle trimmed reads. It also can't handle soft-clipped reads. You might as well not trim. You'll get a lower alignment rate, but your only other choice would be to either exclude the trimmed reads or trim everything down to the same length.
3. You have a mixture of RF and FR alignments and the plotting script can't handle that.
4. Assuming you mean the read length distribution, it's actually shifted to the left, not the right. You had a single length before and a long left tail after trimming.
> 4 votes

---
Tags: rna-seq, rna-splicing, trimming
---


# Question
Title: How to convert a binary matrix of gene presence or absence into a fasta sequence
I have a binary matrix of gene presence or absence which looks like: \[roary output\]

```
Gene    sample1 sample2 sample3 sample4
fliI        1       1       1        1
patB_1      1       1       1        1
pgpA        1       1       1        1
osmB        1       1       1        1
cspA        1       0       1        1

```

How can I convert to fasta so it looks like and technically aligned:

```
>sample1
11111
>sample2
11110

```

Is it possible to do this in R?

# Answer
It's not a fasta file, but:

```
> m
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> # Collapse to labeled strings
> blah = apply(m, 2, function(x) paste(x, collapse=''))
> blah
sample1 sample2 sample3 sample4
"11111" "11110" "11111" "11111"
> # Write that to a file with the appropriate format
> cat(paste(mapply(function(x, y) sprintf(">%s\n%s\n", x, y), names(blah), blah), collapse=""), file="some file.fa")

```

Make sure to change `some file.fa`. I have no idea why you would want to do all of this, but this will give you the output you want:

```
$ cat blah.txt
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

```
> 8 votes

# Answer
A mixture of `apply`, `paste` and `cat` will do this:

```
> data.mat
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> data.conv <- apply(data.mat,2,paste,collapse="") # aggregate along 2nd dim
> cat(paste0(">", names(data.conv), "\n", data.conv), sep="\n")
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

```

If writing to a file is required, just add a `file` parameter to the `cat` command:

```
cat(paste0(">", names(data.conv), "\n", data.conv),
    sep="\n", file="out_fastaLike.txt")

```
> 4 votes

# Answer
One could also use `cut` and a pair of Python scripts:

*transpose.py*

```
#!/usr/bin/env python

import sys

for c in zip(*(l.split() for l in sys.stdin.readlines() if l.strip())):
    sys.stdout.write("%s\n" % ('\t'.join(c)))

```

*mock\_fasta.py*

```
#!/usr/bin/env python

import sys

for e in (l.split() for l in sys.stdin.readlines() if l.strip()):
    sys.stdout.write(">%s\n%s\n" % (e[0], ''.join(e[1:])))

```

Then:

```
$ cut -f2- in.mtx | transpose.py | mock_fasta.py > in.fa

```

Which gives:

```
$ more in.fa
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

```

Why anyone would use R for text parsing is beyond me, but kids these days etc.
> 3 votes

---
Tags: r, fasta, format-conversion
---


# Question
Title: How to convert the .vcf (imputed) file with GT:GP format to GT:DS?
I have the genotyped data from impute2 output in .gen format (imputed to 1000G P3). The file has genotype posterior probabilities (GP:3 values per variant). I have converted .gen to .vcf using qctools and the .vcf file has GT:GP format. I need to convert the .vcf file with GT:GP format to GT:DS. Genotype dosages are recommended for use in qtltools/fastqtl analysis. However, I cannot find a tool that would keep the .vcf format and convert GP to DS. Any help much appreciated!

# Answer
> 5 votes
You can do this in Hail.

Here's the rough code to do it (0.1 versions).

Setup:

```
from hail import *
hc = HailContext()

```

Import the .gen file. VCF works too:

```
dataset = hc.import_gen(
    'src/test/resources/example.gen',
    'src/test/resources/example.sample')

```

Remap the genotype schema and export to VCF:

```
dataset.annotate_genotypes_expr('g = {GT: g.call(), DS: g.dosage()}')\
    .export_vcf('/tmp/out.vcf.bgz')

```

Take a look at the getting started page if you want to try it out!

I should note that you may be able to do QTL analyses in Hail, depending on the method you want to use. See blog post here.

# Answer
> 2 votes
Hm, I didn't know that the plugin existed so I wrote my own script to convert GP to minor allele dosage on github. Maybe someone else will find it useful :) https://github.com/7methylg/VCF-GP-to-DS

# Answer
> 0 votes
There's the dosage plugin for bcftools, but it only outputs tab separated values. It would not be too hard to extend the plugin to output a VCF with the DS tag instead, but it has not been done yet. There's a good chance the bcftools devs would respond to a feature request...

In any case, this code:

```
curl https://raw.githubusercontent.com/samtools/bcftools/develop/test/convert.vcf > convert.vcf
bcftools +dosage convert.vcf > output.tsv
head -2 output.tsv

```

has the output:

```
#[1]CHROM       [2]POS  [3]REF  [4]ALT  [5]NA00001      [6]NA00002      [7]NA00003      [8]NA00004      [9]NA00005      [10]NA00006     [11]NA00007     [12]NA00008    [13]NA00009     [14]NA00010
X       2698560 G       A       0.1     0.0     0.1     0.2     0.3     0.2     0.2     0.2     0.2     0.1

```

This is using bcftools version 1.3.1.

Here is an excerpt from the bcftools manual for the dosage plugin:

> dosage
>
> print genotype dosage. By default the plugin searches for PL, GL and GT, in that order.

---
Tags: file-formats, vcf, format-conversion, impute2
---


# Question
Title: blasting a refseq protein does not show the protein in the result set
Can anyone explain me, why I don't find a specific protein with a blast that was took before from the NCBI refseq database?

Specifically, I was trying to blast the protein with the accession number "NP\_420767" and its sequence, respectively, however that protein does not show up in the results. It not only happens when the standard options are chosen, but also, when "Reference proteins (refseq\_protein)" as database in the blast options is selected.

I am really puzzled by that.. Shouldn't blast show the initial refseq entry in the results list as it is part of the refseq database and has the same sequence?

# Answer
NP\_420767 is represented by the non-redundant refSeq protein WP\_010919826, which has the same amino acid sequence. This is not very clearly annotated, but if you scroll down to the sequence in the GenPept entry for NP\_420767, you'll see the following:

```
CONTIG      join(WP_010919826.1:1..799)

```

See here for more info.
> 2 votes

---
Tags: blast, refseq
---


# Question
Title: Filtering bases based on phred qualities with pysam
Is there a way to filter bases in BAM files based on phred quallities through python's pysam ?

I have a code here that

* Takes the nucleobases per position from a BAM file using pysam's pileup function
* Saves it in ReverseList and ForwardList based on both strands (i.e forward and reverse), I want to reject those bases that have phred quality below 25 to be not stored in the ForwardList and ReverseList lists so that they are not used for further analysis.

```
 samfile = pysam.Samfile( filename, "rb" )
 ReverseList = [''] * lenref
 ForwardList = [''] * lenref
 for pileupcolumn in samfile.pileup() :
     for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
               continue
         if not pileupread.is_del and not pileupread.is_refskip:
               if pileupread.alignment.is_reverse: #negative
                      ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
               else:
                      ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

 samfile.close()

```

Where lenref = 16569(length of mitochondrial genome) and filename is the name of BAM file. I want to filter based of phred qualities of bases.

# Answer
> 4 votes
The PileupRead object as a query\_position attribute, which you can use for this:

```
for pileupcolumn in samfile.pileup() :
    for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
             continue
         if not pileupread.is_del and not pileupread.is_refskip:
             if pileupread.alignment.query_qualities[pileupread.query_position] < 10:
                 # Skip entries with base phred scores < 10
                 continue
             if pileupread.alignment.is_reverse: #negative
                   ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
             else:
                   ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

```

Note the 6th line, which implements a filter with a threshold of 10.

---
Tags: sam, python, pysam
---


# Question
Title: Getting protein FASTA sequence based on keyword with python
I would like to gather proteins FASTA sequence from Entrez with python 2.7. I am looking for any proteins that have the keywords: "terminase" and "large" in their name. So far I got this code:

```
from Bio import Entrez
Entrez.email = "example@example.org"

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)
searchResult = Entrez.read(searchResultHandle)
ids = searchResult["IdList"]

handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
record = handle.read()

out_handle = open('myfasta.fasta', 'w')
out_handle.write(record.rstrip('\n'))

```

However, it can get me several terminases from various organisms, while I need only terminase form bacteriophages (specificly Viruses \[taxid 10239\], host: bacteria). I've managed to get the nuccore accession ids from NCBI of the viruses I am intersted in, but I don't know how to combine those two informations. The id file looks like this:

```
NC_001341
NC_001447
NC_028834
NC_023556
...

```

Do I need to access every gb file of every ID and search for my desired protein in it?

# Answer
> 4 votes
Found what I was looking for. In:

```
searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)

```

I've added:

```
searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)

```

which norrowed down my searches to the desired viruses. Granted it's not filtered by host, but by a taxonomy group, but it is enough for my work.

Thank you @Llopis for additional help

---
Tags: python
---


# Question
Title: Same transcript coordinates in gtf file, different transcript ID
I have a gtf file from Ensembl, and I noticed that several "transcript" annotations have the exact same coordinates. See for instance the third and fourth transcripts ("Y74C9A.2b.1" and "Y74C9A.2b.4") for this gene:

```
$ grep "WBGene00022276" genes.gtf | awk '$3 == "transcript" {print}'
I   ensembl transcript  10413   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.2"; transcript_name "Y74C9A.2a.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS7921";
I   ensembl transcript  11495   16793   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.2"; transcript_name "Y74C9A.2b.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9220";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.1"; transcript_name "Y74C9A.2b.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.4"; transcript_name "Y74C9A.2b.4"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11505   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.3"; transcript_name "Y74C9A.2a.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS16820";
I   ensembl transcript  11618   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.1"; transcript_name "Y74C9A.2a.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS43845";
I   ensembl transcript  11623   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.3"; transcript_name "Y74C9A.2b.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS49800";

```

What does it mean? Do these two "transcripts" actually correspond to an exon, belonging to two different isoforms?

# Answer
`transcript` objects cover the co-ordinates from the start of the first exon to the end of the last exon of a transcript (i.e. an isoform). If two different isoforms share the same first and last exons, but have a different set of internal exons, then their `transcript` entries will be the same, but the set of `exon` entires associated with each transcript will be different.

For example, consider the two transcripts below:

```
                      1         2
             12345678901234567890123
transcript_1 |>>>|----|>>|----|>>>>|

transcript_2 |>>|--|>>|------|>>>>>|

```

These two transcripts would give the following GTF entires

```
chr1   ensembl   transcript    1   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    1   5   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    10   13   .   +   .   gene_id "gene_1"; gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    18   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   transcript    1   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    1   4   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    7   10   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    17   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";

```
> 7 votes

# Answer
As Ian explained, these are different transcripts which happen to have the same start and end positions. You have no information on their exonic structure in that file. However, if you look them up at EnsEMBL, you will see:

Transcript Y74C9A.2b.1 :

```
1   Y74C9A.2b.1.e1  11,499  11,561  -   -   63
    Intron 1-2      11,562  11,617          56
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

```

This has 4 exons, two of which are non-coding (UTR). The second starts at

Transcript Y74C9A.2b.4

```
1   Y74C9A.2b.4.e1  11,499  11,557  -   -   59
    Intron 1-2      11,558  11,617          60
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

```

Granted, the two transcripts are very simlar, they both have 4 exons, the first 2 of which are non coding (UTR). However, the 1st exon of transcript Y74C9A.2b.1 is 63 nt long and ends at 11,561 while the 1st exon of transcript Y74C9A.2b.4 is 59 nts and ends at 11,557.

So there you go: two (slightly) different transcripts.
> 8 votes

---
Tags: sequence-annotation, gtf
---


# Question
Title: Solutions for managing data in a small bioinformatics / 'omics lab?
A different sort of problem: even a small 'omics lab generates a lot of data, raw, intermediate and processed. What (software) solutions exist for managing this data, such that "old" data can be retrieved and checked or re-analysed, even after people have left the lab? Important points would be:

* ease of installation
* ease of putting data in, in an appropriately tagged / labeled fashion (it's no good if the repository is just a centralised bad mess)
* useful search and exploration
* security (i.e. restricted to members of the lab)
* previews / summaries of the data
* can accommodate any dataset
* local, not SaaS

It seems from my explorations of the topic that there's very little between the primitive (e.g. handrolled Access databases) and big industrial solutions. Things I have looked at include:

* Dataverse: very popular, installation seems complex and unclear if uploading is that easy
* DSpace: mostly for publications and documents
* CKAN
* OSF: used this for a while, integrates with a lot of services but uploading data seemed awkward

# Answer
> 4 votes
I don't know of any prebuilt products, but I can describe how we managed this in my postdoc lab, and how I plan to manage it in my newly-started group.

Rule 1: All work happens in the projects project directory on the central filestore, not on your desktop or laptop.

Rule 2: Heavy computational work is done by the groups standard analysis pipelines. Interpretation is done in jupyter notebooks or Rmarkdown.

A project has a directory on the group's filestore. That directory has a fixed structure:

```
proj001----raw_data             *
       |
        --external_datasets
       |
        --src                  *
       |
        --notebooks            *
       |
        --pipeline1
       |
        --pipeline2
       |
        -- etc
 * link to a seperate, backed up filesystem.

```

Pipelines are where heavy duty analysis happens and everyone uses the same standard pipelines, producing the same standard output files and database structure.

So common pipelines might be the mapping pipeline, the readqc pipeline, the differential expression pipeline, the exome pipeline, the chip-seq pipeline etc...

Pipelines have three important outputs: an automatically generated html report, an sqlite database and files in the export directory.

When we archive a project, this is what we save, along with the pipeline's configuration file and log file.

So if I know that Jane Bloggs did an RNAseq 5 years ago on a cell type i'm interested in, if I konw that that was project 5, I know that in the project 5 directory there will be a `diff_expression_pipeline` directory, and that it will contain an sqlite database called `csvdb` and that that will have a table called `refcoding_deseq_gene_diff` and that table will follow a known format. There will be bigwigs in the `export` directory. Or the BAMs will be in the `mapping_pipeline` directory.

Of course the problem remains knowing that Jane Bloggs did this RNAseq and what the project was called. We use a Wiki for this, but its not ideal.

# Answer
> 2 votes
We have been organizing data on a "per project" basis, using a GitHub repository for each project. Each project ends up being a paper (written in R Markdown), so you could have multiple "projects" based on the same datasets. See this as an example:

https://github.com/SchlossLab/new\_project

The analyses are all managed with a Makefile that goes from loading the sequence data from the SRA, to performing analyses, to rendering the final manuscript. This way it is always clear how things were done. Here is an example of what that can look like:

https://github.com/SchlossLab/new\_project/blob/master/Makefile

As far as organizing the raw data and whatnot, I have just used a simple MySQL database hosted on MAMP (Mac version, Linux is LAMP, Windows is WAMP). It has a GUI and everything, and is pretty easy to use. Here is a resource for getting started on that:

http://microbiology.github.io/PDFs/MySQL-Intro-Workshop-Hannigan.pdf

All of this stuff is free to use, which is awesome! Hope this helps! :)

---
Tags: database, data-management
---


# Question
Title: How can I use Nanopore reads to close gaps or resolve repeats in a short-read assembly?
Low coverage MinION reads should be useful to close gaps and resolve repeats left by short-read assemblers. However, I haven't had any success with the software I know about. I'm aware of the following packages, either for scaffolding or closing gaps in short-read assemblies using long reads:

I've tried `npScarf` and `LINKS`, and they ran successfully but didn't resolve the gaps in my assembly. I couldn't get `PBJelly` and `OPERA-LG` to run with the current versions of `BLASR` and `samtools`, and both packages seem not to be maintained. I have not tried `SSPACE-longread` because it's not open source.

**What software can I use to fill gaps and resolve repeats in a short-read assembly with low-coverage Nanopore data?**

---

### More information

I'm finishing a mitochondrial genome. I have a 17 kb short read assembly with one gap. This was made from 2x150b paired-end reads from a TruSeq PCR-free library with a 400 b insert size. Molecular data suggests the genome may be around 32 kb.

After mapping short reads against the short read assembly, `Pilon` identified a tandem repeat of 156 bp in an AT-rich region, but wasn't able to close the gap.

I have 1.4 Gb of MinION rapid reads with an L50 of 8808 b. I mapped these reads against the short read assembly, and I can see reads that span the gap. I seem to have the information required to close the gap, but I don't know how to do it.

The organism has a nuclear genome \> 600 Mb, so long read coverage is low. Despite this, I tried a *de novo* `Canu` assembly of the Nanopore reads, and I pulled out a 39 kb contig containing mitochondrial genes. Most genes on this contig are duplicated and fragmented, and `Pilon` wasn't able to improve it.

Thanks for reading!

# Answer
> 6 votes
You might want to look into Unicycler (manuscript with more information can be found here); even though it is supposed to be used with bacterial genomes only, it might work well with a small genome such as a mitochondrion's.

Beware that if you happen to have very long reads, you might end up with an assembly with multiple copies of the circular genome: you might want to look into circlator then.

---
Tags: nanopore, long-reads, repeat-elements, scaffold
---


# Question
Title: differential gene expression complex design no replicates
We have an experimental design as seen below

Where we administered drug at 0 min for each mouse genotype and took them down at given below intervals. wt and ko mouse models were administered only with the vehicle which would be our negative control.

Performed rnaseq (total rna) truseq stranded protocol would like to see differences between each time point compared to the vehicle, used star aligner to align to mm9 , htseq for gene counts .

```
sample  type    time
WT  Vehicle     30 min
KO  Vehicle     30 min
WT  Drug    30min
KO  Drug    30min
WT  Drug    1 hour
KO  Drug    1 hour
WT  Drug    2 hour
KO  Drug    2hour
WT  Drug    3 hour
KO  Drug    3hour

```

We do not have any replicates and hence cannot use any of the differential expression statistical tools (asked this question at bioconductor where Dr. Love himself replied link to question asked) where he suggested a linear model.

we were wondering if someone could help us with any such design (papers/previously performed analysis) what we would like to see is comparison/differences (significant) between these groups (veh vs 30min, veh vs 1 hour, beh vs 2 hours, veh vs 3 hours)

# Answer
> 7 votes
Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.

# Answer
> 7 votes
Mike Love is right. **If** the response you are looking for is linear in terms of change per minute, the most productive approach is likely to be fitting a linear model. You **might** get something out of this because the difference between successive time points represents multiple measurements of the rate of change over time. The biggest problem is that the time points aren't independent. You would need to fit a model that estimated that rate for both WT and KO and then test the difference in rates.

I see three ways forward, in all cases you want to test the last term in the model `~time + genotype + time:genotype`

1. Ideally you would like to fit a negative binomial linear model to the read counts. But that is going to be difficult as you have no replicates from which to assess the dispersion. Thus you could generate read counts using your favorite read count pipeline (STAR+featureCounts/Salmon/Kallisto) and take the approach in the edgeR guide noted by b.nota and then use edgeR to test the fitted model.
2. You could use something like salmon or kallisto to estimate the TPMs assume that log TPMs are normally distributed and model them with plain `lm` in R, using something like `broom` to model each gene separately. A tutorial by David Robinson, the author of `broom` on doing this sort of thing can be found here
3. You could use the Salmon or Kallisto bootstraps to estimate *technical* variation (but not biological) and then I suspect you could use Sleuth to fit the model.

In any case these results should be treated with scepticism andmostly treated as hypothesis generating. You might find genes that you want to go back and confirm and you might see some patterns of interest if you look at some enrichments etc. which you can then design experiments to test.

---
Tags: rna-seq, deseq2, differential-expression
---


# Question
Title: How to safely and efficiently convert subset of bam to fastq?
### Question

How can I extract reads from a `bam` file (produced by `bwa-mem`) to `fastq` given a list of reference sequences to filter out?

### Potential difficulties

* maintaining FR orientation of pair end reads (in `bam` all the sequences are reference sequences)
* keeping R1 and R2 reads
* keeping quality scores in the same encoding as original fastq (default illumina phred scores in my case)
* bam can be (ana usually is) sorted by coordinates

### Almost there solutions

1. Sujai's perl solution in blobology does exact opposite - getting reads from list of references (so I could just reverse the list). The disadvantage is that the script outputs an interleaved `fq` file; requires unique names of mates, otherwise R1/R2 information is lost.
2. samtools + grep them all from fastq files

create a list of read names that do not map to filtered scaffolds. (cut will extract just read names, uniq will collapse pair end read names if they are the same). Then grep read names from fastq files and remove -- separator between hits

```
samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | cut -f 1 | uniq > list_of_reads_to_keep

grep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^--$" > foo_R1_filtered_bash.fq
grep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^--$" > foo_R1_filtered_bash.fq

```
3. filter bam & picard tools

Or I could do just the filtering part and use Picard-tools (Picard.SamToFastq), but as usual I am avoiding java as much as I can. I guess

```
samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | java -jar picard.jar SamToFastq INPUT=/dev/stdin \
  FASTQ=foo_R1_filtered_bash.fq SECOND_END_FASTQ=foo_R2_filtered_bash.fq

```

---

The first solution does not really work for me since I do not want to rename all the reads in bam file and I want to keep the R1/R2 information (since R1 and R2 have different error profiles). Both solutions 2 and 3 I find bit clumsy and I am not sure if they are general, I might get some unexpected behaviour if one of reads is mapped second is not... They both relay on the same filtering step.

I was wondering about some pysam solution. I guess it will be much slower, but at least it will be much clearer and perhaps more general. Something like in Convert Bam File To Fasta File \- there is pysam solution for fasta (not fastq), almost there...

### Story

I have very fragmented reference genome. Some of scaffolds them are too small to works with, and some of them are contaminants (identified using blobtools). I want to separate reads that are mapping to different groups to separate contaminants, short scaffolds and scaffold that will be used for downstream analysis. The reason is that if we remap all the reads to filtered reference (0.7 - 0.8 of original genome), the most of them (0.95 - 0.99) will still find a place where they map, therefore there is 0.2 - 0.3 of misplaced reads that will obviously have to affect downstream analysis, like variant calling.

This filtering idea is based logic that if the filtered duplicated genomic region will contain some small differences, they will attract their reads (and if I filter them I will improve variant calling) and if they will be exactly same, they will get reads assigned at random, so there is no harm in doing that.

# Answer
I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter list), store the properly oriented sequence/quality/read name in a buffer, and then dumps that buffer entry to disk once the mate is found. This should be reasonably performant (hey, it's python, to don't expect too much), though if you happen to have indexed BAM files then one could think of ways to make this run faster.

Do check the output, since I've only run one test.
> 5 votes

# Answer
Ok, I wrote a bit brute pysam / BioPython parser that uses index of bam to get proper order of read pair for R1 / R2 files and bitwise flag. It should not be too difficult to add more sophisticated filtering rules now.

```
#!/usr/bin/env python3
# 1. arg - indexed bam file
# 2. arg - list of headers to filter
# 3. arg - name pattern for the output reads

import os
import sys
import pysam
from Bio import SeqIO, Seq, SeqRecord

samfile = pysam.AlignmentFile(sys.argv[1], "rb")
header_set = set(line.strip() for line in open(sys.argv[2]))
base = sys.argv[3]

out_R1 = base + 'R1_filtered.fq'
out_R2 = base + 'R2_filtered.fq'

with open(out_R1, mode='w') as R1, open(out_R2, mode='w') as R2:
    for entry in samfile:
        if entry.is_read1 and not entry.reference_name in header_set:
            # pait pair
            entry_R2 = samfile.mate(entry)
            # do some sequence gymnastics with R1
            seq_R1 = Seq.Seq(entry.seq)
            if entry.is_reverse :
                seq_R1 = seq_R1.reverse_complement()
            # do some sequence gymnastics with R2
            seq_R2 = Seq.Seq(entry_R2.seq)
            if entry_R2.is_reverse :
                seq_R2 = seq_R2.reverse_complement()
            R1.write('@' + entry.qname + '\n' + str(seq_R1) + '\n+\n' + entry.qqual + '\n')
            R2.write('@' + entry_R2.qname + '\n' + str(seq_R2) + '\n+\n' + entry_R2.qqual + '\n')

samfile.close()

```

There are some ugly bits, please feel free to make it nicer.
> 2 votes

---
Tags: sam, fastq, format-conversion, filtering
---


title: Fasta Sequence Identifier format?

text:
# Question
Title: Fasta Sequence Identifier format?
I have been analyzing some virus DNA from the NCBI databases. The fasta sequences that I recieve from them have header lines that look like this:

```
>gi|61393989|gb|AY848686.1|

```

I think that the second number is the GenBank ID, but I have several questions:

1) What does the "gi" in the first ID stand for?

2) Is this a standard format for storing a series of IDs associated with sequences? That is, should I expect to see this format elsewhere and be prepared to parse fasta headers like this in code that I write?

3) It seems that `AY848686` is also a valid GenBank ID, and I see the `.1` missing/present in GenBank IDs in various places. What does the `.1` mean at the end of GenBank IDs and why is it there/does it *need* to be there?

Thanks in advance for you answers.

# Answer
> 5 votes
1. The `gi` is an abbreviation for "Genbank identifier".
2. This is a pretty standard convention used by data stored in NCBI databases. There used to be a pretty comprehensive description of the conventions used at NCBI (I wouldn't say it was a standard or specification, just convention) here, but this page is no longer available it seems.
3. The `AY848686` label is an accession value, and the appended number is a version number.

See this page for more info.

The short answer is: these are all internal conventions. None of it *needs* to be there, except perhaps to ensure NCBI's internal tools work as expected.

---

**UPDATE**: I found a non-official copy of the old docs here. Look for the section starting with "FASTA Defline Format". I've copied it here to prevent further problems with bitrot in the future.

```
FASTA Defline Format
--------------------
The syntax of FASTA Deflines used by the NCBI BLAST server depends on
the database from which each sequence was obtained.  The table below lists
the identifiers for the databases from which the sequences were derived.

  Database Name                         Identifier Syntax

  GenBank                               gb|accession|locus
  EMBL Data Library                     emb|accession|locus
  DDBJ, DNA Database of Japan           dbj|accession|locus
  NBRF PIR                              pir||entry
  Protein Research Foundation           prf||name
  SWISS-PROT                            sp|accession|entry name
  Brookhaven Protein Data Bank          pdb|entry|chain
  Patents                               pat|country|number
  GenInfo Backbone Id                   bbs|number
  General database identifier           gnl|database|identifier
  NCBI Reference Sequence               ref|accession|locus
  Local Sequence identifier             lcl|identifier

For example, an identifier might be "gb|M73307|AGMA13GT", where the "gb" tag
indicates that the identifier refers to a GenBank sequence, "M73307" is its
GenBank ACCESSION, and "AGMA13GT" is the GenBank LOCUS.

"gi" identifiers are being assigned by NCBI for all sequences contained
within NCBI's sequence databases.  The 'gi' identifier provides a uniform
and stable naming convention whereby a specific sequence is assigned
its unique gi identifier.  If a nucleotide or protein sequence changes,
however, a new gi identifier is assigned, even if the accession number
of the record remains unchanged. Thus gi identifiers provide a mechanism
for identifying the exact sequence that was used or retrieved in a
given search.

For searches of the nr protein database where the sequences are derived
from conceptual translations of sequences from the nucleotide databases
the following syntax is used:

                     gi|gi_identifier

An example would be:

        gi|451623           (U04987) env [Simian immunodeficiency..."

where '451623' is the gi identifier and the 'U04987' is the accession
number of the nucleotide sequence from which it was derived.

Users are encouraged to use the '-I' option for Blast output which will
produce a header line with the gi identifier concatenated with the database
identifier of the database from which it was derived, for example, from a
nucleotide database:

        gi|176485|gb|M73307|AGMA13GT

And similarly for protein databases:

        gi|129295|sp|P01013|OVAX_CHICK

The gnl ('general') identifier allows databases not on the above list to
be identified with the same syntax.  An example here is the PID identifier:

        gnl|PID|e1632

PID stands for Protein-ID; the "e" (in e1632) indicates that this ID
was issued by EMBL.  As mentioned above, use of the "-I" option
produces the NCBI gi (in addition to the PID) which users can also
retrieve on.

The bar ("|") separates different fields as listed in the above table.
In some cases a field is left empty, even though the original
specification called for including this field.  To make
these identifiers backwards-compatiable for older parsers
the empty field is denoted by an additional bar ("||").

```

# Answer
> 1 votes
1. That is a GenBank gene ID and it (apparently) stands for "GenInfo Identifier". The full GenBank ID is `gi|61393989`.
2. Yes, this is the standard format for sequences provided by NCBI. However, you can have all sorts of extra details appended to that. They (NCBI) generally use `|` as a field separator, so you could try to write code that parse that, but I really wouldn't recommend it. The FASTA headers have traditionally been free. People can put anything they like there apart from a `\n`, so I really wouldn't put any effort into writing code to parse these headers.

If you need a short name for the sequence, just take the first N characters or, if it matches `^..|.+?|`, take the characters until the second `|`. Seriously though, don't bother parsing this it is very rarely helpful and is almost certain to break as soon as you use a different database.
3. The AY848686.1 is a GenBank accession. The `gb` indicates that it's a GenBank accession and the `.1` means it is the 1st version of that accession.

As for why there are two (source):

> Unlike the gi number system, in which sequence identification numbers were not necessarily consistent across the databases (e.g., GenBank and EMBL could each assign their own gi number to a sequence), the new system is designed to ensure consistency. It is also designed to show a relationship between a sequence identification number and the accession number of the record in which it is found. In contrast, GI numbers are assigned consecutively and bear no resemblance to the accession number. Finally, the new system allows the assignment of alphanumeric protein IDs to proteins translations within nucleotide sequence records. The protein IDs contain three letters followed by five digits, a period, and a version number.

---
Tags: fasta, identifiers
---


# Question
Title: Is there a standard k-mer count file format?
I am doing a research project involving calculating k-mer frequencies and I am wondering if there is any standard file format for storing k-mer counts.

# Answer
> 6 votes
Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the `dump` command (but its default output format is a custom binary format):

> The dump subcommand outputs a list of all the k-mers in the file associated with their count. By default, the output is in FASTA format, where the header line contains the count of the k-mer and the sequence part is the sequence of the k-mer. This format has the advantage that the output contains the sequence of k-mers and can be directly fed into another program expecting the very common FASTA format. A more convenient column format (for human beings) is selected with the -c switch.

Jellyfish changed their internal format between v1 and v2 (both not FASTA), because they changed to doing counts based on bloom filters. Jellyfish2 has an optional two-pass method that sets up a bloom filter intermediate file to record kmers, and multiple different final reporting formats.

Khmer also uses bloom filters, but in a slightly different way. It also has been extended to be useful for partitioning and comparing datasets.

# Answer
> 5 votes
Based on my experience reviewing k-mer counting software and as a core contributor to the khmer project, I can confidently say that there is no widely used standard format.

---
Tags: file-formats, k-mer, phylogenetics
---


# Question
Title: When performing differential expression analysis, should genes with low read counts be removed before or after normalization?
I have RNA seq data which I've quantified using Kallisto. I'd like to use tximport to transform the read count data into input for EdgeR, following the R code supplied in the tximport documentation:

```
cts <- txi$counts
normMat <- txi$length
normMat <- normMat/exp(rowMeans(log(normMat)))
library(edgeR)
o <- log(calcNormFactors(cts/normMat)) + log(colSums(cts/normMat))
y <- DGEList(cts)
y$offset <- t(t(log(normMat)) + o)
# y is now ready for estimate dispersion functions see edgeR User's Guide

```

I would ideally follow this up with EdgeR's exactTest. However, I'd also like to remove genes with low read counts using code such as this:

```
filterCounts <- function(counts) {
  cpms <- cpm(counts)
  keep <- rowSums(cpms > 1) >= 1
  counts[keep,]
}

```

The only issue is that I can't tell whether I should remove the genes with low read counts before or after I normality the counts while importing with tximport in the first code. I'm thinking that it would be best to do this before, editing out the abundance, counts, and length rows of the txi data frame which correspond to lowly expressed genes, but I would like this to be confirmed before I proceed.

# Answer
> 6 votes
The more genes you have the more robust the scaling factor is (among the reason why one doesn't normalize to ERCC spike-ins without a compelling reason), so I suppose in theory it's better to filter after determining the scale factor. Having said that, I'd be surprised if the results changed much either way. Unless you end up filtering out a LOT of genes the scaling factors should be fairly robust to the change.

---
Tags: rna-seq, statistics, edger, best-practice
---