text:
How to filter info score post-imputation?
this is a very basic question, but I cannot find it explicitly stated anywhere: When exactly should I filter the info scores after imputation?
I did imputation in Impute2 and am planning on using SNPtest to analyze, but cannot tell if I should filter the scores in Impute2, SNPtest, Gtool, QCtool, or through python.
I'm hoping to get rid of the poorly imputed results (with an info score cutoff <0.5) in a population of ~4000 (half cases, half controls). After this, I'm running an association analysis -- probably through PLINK, but I'm still figuring out the best way to go from Impute2 to PLINK.
Any help would be greatly appreciated!!

In simple GWAS setups, each SNP is analyzed independently.  In those cases you can filter out SNPs with poor INFO scores at any point.
For analyses that combine information across SNPs (for example FaSTLMM), I would recommend filtering out SNPs before running the step that combines information across SNPs.
So, if you filter right after running Impute2, you should be safe.

Thanks so much, winni2k. What is the best way to filter the SNPs? Is that usually an option when inputing files to a software, or something that should be done separately (e.g., through something like python)?

Plink is wildly popular and very fast. But do use version 2 if you want to avoid surprise allele reordering.

Depending on the size of your samples (1K, 10K, 100K?) or the size of imputation region (1 gene, 1 chromosome, whole genome?) imputation files even after zipping are usually huge.
Instead of filtering SNPs and creating yet more gazillions of subset files. I prefer to have excludeSNP.txt file list of SNPs. This list then used as input for further analysis using SNPTEST flag -exclude_snps_g. This list is part of IMPUTE2 output or could be additional list of SNPs that we wish to exclude for other reasons. In short, filter at the point of analysis not the imputated files.
If info file is missing we can run SNPTEST with -summary_stats_only flag, which gives you the info score.
Here is an old still relevant post at BioStars post-imputation QC.


text:
Basecalling process carried out by full_1dsq_basecaller.py/.exe
I am unable to understand what exactly the full_1dsq_basecaller.py/.exe script does. Is it giving the consensus sequence for linear and complementary strand, or it is just giving the reads that came from the complimentary strand?
If we are getting two fastq files as output which one we should use for alignment?
Thank you in advance.
I greatly appreciate any help

Possible duplicate of Removing repeated reads from nanopore 1D² reads

My understanding is that the 1D$^2$ process is carried out in three stages:

Carry out a 1D calling of reads (i.e. what the standard Albacore does at the moment)
Link reads that are similar to each other (in reverse-complement orientation), from the same channel, and have a very short pore transition time
Create a 2D consensus base call, using information from the two linked reads

The two outputs from the basecaller are the 1D calls from step 1, and the consensus calls from step 3. Ideally just the 2D consensus calls should be used for assembly (assuming it represents a random distribution of reads), but  if that gives too low coverage (or is not random enough), then the non-converted reads from the 1D calls could be included as well.
Because the linking is not a physical link, it's not a good idea to use 1D$^2$ for amplicon reads; the likelihood of two complementary reads from different templates appearing one after another is too high.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
How to summarize multiple exon copy numbers into copy number of the corresponding gene
I have a matrix, sample by exon, containing a copy number value for each pair (sample, exon). I would like to generate a second matrix, sample by gene, where the copy number of the exons is consolidated in a single value per gene.
Due to the nature of copy number alteration, I am not sure if this is meaningful, or feasible (e.g., how can it be summarized, in the case a gene is interested by 2 alterations with different copy numbers?).
Is there a common approach that is used in this case?
Edit
A possible example is the following:
         gene1_exon1  gene1_exon2  gene1_exon3 ... geneN_exonM
sample1            2            1            1               0
sample2            2            2            1               0
...
sampleJ            3            2            2               2

A possible summarized matrix (with majority voting - not sure if this would be the best approach):
         gene1 ... geneN
sample1      1         0
sample2      2         0
...
sampleJ      2         2

I believe this is more an issue related to the biological interpretation than the actual implementation (i.e., I am already implementing a function in pandas). My concern is on the interpretation, since my prototype select as gene copy number the majority copy number from its exons. However, this has some caveats (e.g., what if you have a GAIN on an exon, and a LOSS on a second exon, for the same gene - will you call GAIN or LOSS?). The big picture is that I want to compare two samples CNAs at the gene level, but I could not find this task discussed online.
The main objective of this summarization is to compute gene level concordance between samples. Exon-level or base-level concordance is computed in parallel, in order to have different CNA granularities.
However, one might not use the majority voting proposed and, e.g., call a CNA on a gene if the alteration is found in at least p% of exons (e.g., 30% of exons) - this is just an example, I am not suggesting to use this approach.

You might be interested in the approach we took when trying to find CNVs that overlap known, reported ones. have a look at our documentation here: https://varsome.com/about/resources/sv-implementation/. I'll try and get the relevant bits into an answer if I get the chance.

Thank you! I will look into it

Take it with a pinch of salt, of course: that was not peer reviewed and was just the best system we could come up with based on our own best judgment. That said, I do think it is a reasonable approach, I am just aware that I have been wrong before! :)

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
What are RIKEN genes?
I am using data from the Tabula Muris Consortium.
These are gene counts from scRNA-seq of mouse cells.
There are some specific genes with name ending with suffix 'Rik' (e.g. 0610005C13Rik or 0610007C21Rik). Searching online I found that they are from RIKEN project.
However, they are different in some respects (e.g. antisense lncRNA gene, or mouse 0610007C21RIK Adenovirus).
I want to exclude them from downstream analysis, but it is unclear to me the main purpose of this collection. For which kind of analysis are they used?

My understanding is that these are cDNA sequences from the RIKEN project that don't yet have another gene name. The project carried out a comprehensive sequencing of the mouse genome, and as a result found a lot of cDNA transcripts that were previously unknown. Now they have some annotation (which can be varying degrees of functional annotation), and an accepted official identifier to allow people to talk about these transcripts when discovered in other studies.
It'd be interesting to see if any novel transcripts have been found from the Tabula Muris project. They would also need their own identifier, possibly leading to similar intrigue about their history.

cDNA is generated from RNA so they most probably sequenced RNA (the transcriptome) instead of the genome itself.

To extend a bit more to @gringer answer. RIKEN is a large Japanese research institute, that participates in different scientific collaborations.
One of these is the FANTOM collaboration (Functional ANnotation Of the Mammalian genome), which lately have resulted in the release of FANTOM CAT, a browsing tool for exploring lncRNA relative to associated Genes, Traits and Ontologies.
Many of the transcripts and genes discovered during this work, are not yet named, and are therefore called CAT-genes, following numbering comparable to ensembl (eg. CATG00000118335.1).
I'm unsure how RIK genes have been discovered, as I have not personally encountered these genes yet, but in case you likewise encounter more genes discovered in collaboration with RIKEN, you now also know of CAT genes.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Web scraping in R
I have ~130 gene names of different human proteins.
I'm looking for a convenient and systematic way to locate each gene's location in the genome, extract its CDS sequence and find whether there are SNPs in the CDS (in any context) and what's their position.
Moreover, I've chosen for each of those proteins, several amino acids for whom I'm trying to find if the SNP fall in their exact coding nucleotide region (the amino acid triad of nucleotides).
Alternatively, I'm looking for an algorithm that would take the gene name out of my list, find its location in the genome and look for SNPs in the exons - if it would found any SNPs, then it should translate the exon to an amino acid sequence and find whether one of my pre-defined amino acids from each protein is translated by the nucleotide site (triad of nucleotides) with the SNP.
I'm using R, and I thought to use web scraping (should I use the rvest package?) but I'm not quite familiar with the technique. Is it possible?
What websites should I use? Genome Browser? EMBOSS Transeq? Other websites?
Thank you

Why do you need to do web scraping instead of using biomaRt and the packages of Bioconductor to find and retrieve the information? Did you look up how do other papers obtain the sequence and positions? I think your attempt parsing webs is full of problems if not unfeasible.

"convenient and systematic way to map genes..." we can only dream.

You can search databases like Cosmic(Cancer) with the gene names to get frequently occurring mutations in cancer. I'm not sure about any specific tool/pipeline available for this exact scenario.
As an alternative approach literature mining for gene associated mutation discovery, you can use text mining with machine learning to automate the whole process.
With no specific context in mind, you will end with several combinations of mutation per gene. Then for each variation translate the DNA sequence and match the result with reference peptide sequence you have.

I don't have a complete answer to your question, but as a primer I can tell you that you should look into the Bioconductor packages. The reason is that biological sequence data are deposited in online databases (like GenBank, and others), and exposed by these databases in such a way that you precisely don't need web scrapping techniques to retrieve the data you're interested in.
rvest is a great package for web scrapping, but definitely not the appropriate tool to solve your problem.


text:
Is there a tool that can take a protein's amino acid sequence and would display it's locus on the genome?
I have the UNIPROT IDs, PDB IDs and FASTA files of several known proteins. I am looking for a tool that can take as input the protein's amino acid sequence and display the coding nucleotides of those amino acids. Moreover, I look for a tool that would be able to display those nucleotides' location on the human genome.
The final outcome I desire to achieve is to find SNPs in the exons of the protein.
A tool that can handle several amino acid sequences at once is preferred for I have plenty of proteins. I've tried to use UCSC's BLAT human genome search but it was too cumbersome. I've also tried BLAST, but I'd like a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.
Thank you for your assistance.

For something "straightforward" you'd need the correct reading frame, strand orientation, and splicing to get back to a nucleotide sequence to align to the genome. The tools mentioned should help take a lot of work out of this.

I think the mistake here is thinking that you need to use the sequence in order to find the genome co-ordinates. If these proteins come from a well annotated speices, you can easily do this with BioMart.

Go the the ENSEMBL web site.
From the top menu, select BioMart
Under "-CHOSE DATABASE-" select Ensembl Genes, and under "-CHOOSE DATASET-" select your species (eg. Human Genes).
Select "Filters" from the lefthand menu and then expand the "GENE" section.
Select the "Input external references ID list [Max 500 advised]" and from the drop-down select "PDB IDs" or the correct "Uniprot KB" ID.
Paste your IDs into the box.
Select "Attributes" on the left-hand menu. Making sure you've the "Fetaures" radio box select, you can select Chromosome and either Gene Start/Gene End or Transcript start/Transcript end from the GENE section. These coordinates will also contain the untranslated regions unfortunately.  You may also want to play about with the various things under "PROTEIN DOMAINS AND FAMILIES", particurlarly under "Protein Features", although I'm not sure what they all mean.
When you are happy with the selection, select the big "Results" at the top to get the coords.

You can also get out Sequences for both the full transcript and the CDS and known germline and somatic variants that fall within the boundaries of the genes.

Thank you, very useful tool. Can it display the external results references IDs (ie the PDB IDs)  at the results page?

Have you tried BLAST, or specifically tBLASTn?  This will take a protein and compare it any set of reference genomes.  You can also BLAST multiple sequences at once.

I've tried it, but I look for a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.

I think that you will have to define "straightforward" if you want help.  BLAST and BLAT are the standard tools.  You need to describe what issues you have with those tools and what you would like instead.


text:
Test for differences between groups of samples
Sorry if the answer to this should be obvious.
I have RNA-expression results from 24 samples which can be divided into 6 groups, (wildtype and two different mutants at two different ages) with a total of 24 samples. I have made a PCA plot and when describing it in my thesis said separation between samples is minimal. One of my examiners has said

You state that the separation between samples in the principal component analysis is minimal, this needs to be quantified to differentiate between minimal but statistically significant, and not significantly different.

I have no idea how he expects me to test for significance. I know how to test for significant differences between individual samples, and for differences per gene/principle component. I don't know how to group samples together while testing for significance against multiple things or if I should use genes or principle components (probably the latter).
I'm writing it in R so if anyone knows which R function can do the kinds of statistical tests I need that would be a huge help. Alternatively is this something which I shouldn't really do and I should push back on?
I'm already performing significance testing using DESeq2. To the best of my knowledge that only tests for whether a single gene is significant. It doesn't test for whether the samples are significantly different as a whole which is what I'm after.

Hi Sethzard I think you need to do is present the experimental design. The problem at present is it's not clear what the criticisms are because we don't know your experiment for "grouping samples together". @RamRS is correct you cannot use PCA. At least you cannot use traditional PCA as used in RNAseq. It's also clear "I don't know how to" is this code or test?

PCA does not do any testing and has no p-values. You may state that a separation of less than m units in PC > i is not significant while a separation of greater than n units in PC < i is significant as PCs < i explain X% of the variation (i being 2 or 3). That's how you can back your statement on the separation being minimal.
That is, in this sample plot:

PC2 distance of < 0.05 might be called trivial compared to PC1 distance > 0.5 if PC1 explains 80% of the variance. This is just 2 PCs though, so your actual case will need to be more robust.

I know that PCAs don't do any testing on their own. I wasn't sure if they could be used as part of a statistical test which does test for significance or if I'd need to use the raw values. Thank you for the suggestion of other ways to do it.


text:
Large range parameter identification
I have a model (ODE) of an enzymatic reaction, and I'd like to estimate four of its parameters to fit some available data. To obtain the model I am just using Mass-Action Kinetics.
Some more details:

The model is of the form $S_1+2S_2+S_3+E\leftrightarrow C \leftrightarrow P_1+P_2+P_3+E$, where $S_i$ are substrates concentrations, $P_i$ products concentrations, $E$ enzyme concentration, and $C$ the complex concentration.
$S_2$ is Hydrogen and $P_3$ is Ammonium ($NH_4$), so besides the aforementioned reaction I have one reaction for the buffer and one reaction for ammonium dissociation, however, the parameters involved are known.
Substrate concentrations are in the range of $1-10mM$.
Enzyme concentrations are in the range of $.1-10\mu M$.
The parameters I want to estimate are the 4 reaction rates, say, $k_1,k_{-1},k_2,k_{-2}$.
The available data (measurements) is only pH. Thus I have a simulation of the model where the final output is $pH=log_{10}(1/S_2)$.
Then, I try to estimate the parameters from the data and the (simulated) pH.

To estimate the parameters I am using Matlab, more specifically the 'Identification Toolbox' and the 'lsqnonlin' function.
I set the initial guess "visually", that is, I move manually the parameters and whenever I assume the simulation is "close enough", I set the initial parameters guess.
The problem is that some of the parameter of he initial guess are of order $1E20$. I am thinking that because of the large value of the parameters, the algorithm cannot find any better parameter values. I have tried moving several optimisation options without any valuable result.
Q: Is there some specific software, or method that is commonly used to estimate parameters of biochemical reactions? Is the large values of kinetic rates a common problem in parameter identification for these type of systems?
I do not know any $P_i$ though, just pH as a function of time. In principle nor E neither C are constant, but I could assume quasi-steady state conditions, no problem. My feeling is that I am not using or setting up correctly the numerical tool...


text:
How to compare transcriptomic profiles of two cell types (single cell RNA-seq)?
I found this interesting Single RNA-seq data set in GEO, but I am not sure how to analyze it appropriately.
They have deposited transcriptomic profiles of human and mouse pancreatic islets (pancreatic cells: Beta cells, Delta, etc).  The problem I see is that the different pancreatic cell types are not in equal numbers: 872 Beta cells and 214 Delta cells.
I want to see whether the expression of a set of genes are correlated between the two cell types. I don't think that I can visualize the correlation between the expression of the Beta and Delta cells on a scatter plot if they are unequal in number. Is there a workaround to include all the cells in the analysis, without losing good quality transcriptome of beta/delta cell?
What I would like to do is to compare the expression of two pancreatic cell types (Beta v.s. Delta cells) using a scatter plot. I have a set of genes that I am interested in, and would like to see their expression profile for beta and delta cells. As far as I know, it is only possible to visualize this using a scatterplot (X:Beta vs Y:Delta) where the number of isolated beta cells is the same as the number of delta cells.
Question: Given unequal number of isolated panacriatic cells, what would be an appropriate way to compare the expression profiles? Should I just ignore the extra ones?
Any idea how to approch this problem?
Thanks,

What is wrong with a traditional differential expression analysis? How big is the difference in number between Beta and Delta cells?

@Llopis Quite large ... the difference between beta and delta is 658 cells.

Do you have a particular set of genes in mind or are you hoping to go fishing for one? It's also unclear why you can't do visualizations with unequal sample sizes.

Any clarifications (such as what you have just written) should be put in the question, rather than the comment.

What makes you think you wouldn't be able to visualise the correlation for unequal numbers?

@MEhsan You can randomly pick the same amount of delta cells than of beta cells and bootstrap to find the real correlation (or plot the scatter plot)

@Llopis I'm not familiar of how to use bootstrapping and randomly sample the cells.  Could you give me an example?  Is there a way to do this in Python ...?

It is "just" a for loop where you pick different samples, I have never done it in Python so you'll need to search in the web.

You are looking for a way to compare expression profiles between cell types, and your data is counts of genes expressed per cell. Your issue is that the data is highly dimensional. It has as many dimensions as you have genes. To be able find useful information, you must reduce these dimensions using a statistical technique such as Principal COmponents Analysis (PCA), t-Distributed Stochastic Neighbor Embedding (t-SNE), or Multidimensional Scaling (MDS). The result of these analyses is a reduction of your thousands of dimensions to vectors (eigenvectors in PCA) of variance which you can plot against each other in a scatter plot.
I have been working with Single Cell RNA-seq data and have been exploring the difference in gene expression in a certain type of cells between timepoints. Seurat (R package) is a very powerful tool for scRNA-Seq data analysis.
Here is a a tSNE of my data. I have added the timepoint labels, but you can see that after performing the dimensionality reduction the two different timepoints cluster by their gene expression profiles. I can then use Seurat to see which specific genes are differentially expressed between the clusters and a whole bunch of other fun things.

I just add the module to run t-SNE in python, since OP seems to want to use a python module: http://scikit-learn.org/0.15/modules/generated/sklearn.manifold.TSNE.html

Scatterplots are used to represent multiple dimensions of the same piece of data. A scatterplot would be used in this case to display all genes in your set on a single plot, but I would have some difficulty interpreting such a plot. If you want to look separately at each gene, then a density curve (or histogram) would be more appropriate.
Your problem is similar to common problems in flow cytometry, and I don't think there's any reason why similar approaches can't be used. Here's a histogram / density plot showing two different types of cells (negative control, and treatment condition), indicating the difference in fluorescence of a particular marker (somewhat analogous to gene expression) across all cells of each type. Note that the number of cells are different between the negative control and treatment:

This image was taken from here, which also demonstrates a 2D scatter plot (where each point is a cell, rather than a gene).

seurat3 may help. You can compare the marker genes from the 2 cell types or treat one cell type as a control. Please see the example analysis from: https://satijalab.org/seurat/v3.0/immune_alignment.html


text:
DeSeqDataSet experimental design: column with integer values
I am creating DeSeqDataSet and I am unsure whether I need to include the number of mice into the design formula. Basically, the samples look like this:

I have 4 different groups: GF, B6, ET and Wld. It is obvious that these should go into the design formula, but I am not sure about the mice number. Basically for each group, e.g. GF, there are 4 experiments with 1, 2, 3 and 4 mice each. When I am trying to use both Groups and Mouse columns as a design formula, R is telling that I need to convert Mouse to factor:

I am a bit confused about two things here: first, as I told already, should I include this Mouse column at all in the design formula of DeSeqDataSet and second, if I do need it, how to convert it to factor? Any suggestions would be greatly appreciated.
Update
It turned out that I was wrong in my initial assumptions: Mouse# really was an ‘id’ which I did not include in the design formula. There were 4 groups of 16 mice in each for 4 trials.
What concerns me also is that if I specify design formula as ~ Groups + Mouse, the DeSeq() function starts comparing the Mouse but not the Groups which is the most interesting for me.

Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples.
In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. This would also be the case if the same mouse were being used for multiple RNASeq experiments. Converting an integer value to a factor is done by the factor method:
sampleTable$MiceUsed <- factor(sampleTable$MiceUsed);

Other examples of things that might be useful in the design formula would be sex, age, and sequencing run batch (assuming these are not consistent across samples).
With a design formula of something like ~Groups + Mouse, it will also be comparing Groups, it's just hidden in the results. The DESeq2 vignette is really great at explaining what's going on, but as a quick answer you can pull out the Groups results with the contrast function. This will pull out group differences that are independent of mouse number:
results(dds, contrast=c("Groups", "GF", "Wld"))

[first value is the variable name, second value is the numerator condition, third is the denominator condition]


text:
Large dataset normalization for PCA
I need to normalize a large (400Mb) dataset for doing PCA analysis. I want to use scran for doing that:
biocLite("SingleCellExperiment")
biocLite("scran")
library(SingleCellExperiment)
library(scran)
list_of_sce <- list()
# Looping though the UMI_count 'split_factor' columns at a time
split_factor = 500
for(i in seq(1,ncols, split_factor)) {
  num_loop = floor(i / split_factor) + 1
  idx = ncols
  if (i + split_factor < ncols) {
  idx = i + split_factor
  }
  sce <- SingleCellExperiment(list(counts=UMI_count[, i : idx]))
  # Normalization of the dataset containing
  # heterogenous cell data (different cell types)
  clusters <- quickCluster(sce) ## <- ** error appears here **
  sce <- computeSumFactors(sce, cluster=clusters)
  list_of_sce[[num_loop]] <- sce
}
cbind(list_of_sce)

However, when I try to run this it produces an error at the quickCluster step:
Error: cannot allocate vector of size 23.7 Gb

It feels like cbind at the end is not the right way to go, and I need to somehow combine sce and re-normalize them again somehow. The question is how to do that in a correct way: normalizing large dataset for PCA analysis when it can not be fully loaded into memory?

BTW, there are many datasets on kaggle.com

Which error do you get when you normalize the whole dataset? How big is your dataset?

I normalized a high throughput dataset for a school project using DESeq library using the script bellow. The code is based on a lesson I had. My goal was determine the over expressed genes, but the normalization step should be the same.
#####################################################################
### DATA PRE-PROCESSING: Normalising gene expression distributions
#####################################################################
# Genes that are not expressed at a biologically meaningful level
# in any condition should be discarded to reduce the subset of genes
# to those that are of interest, and to reduce the number of tests carried out
# downstream when looking at differential expression.
# Using a nominal CPM value of 1 (which is equivalent to a log-CPM value of 0
# a gene is deemed to be expressed in a given sample if its transformed count
# is above this threshold, and unexpressed otherwise.
# Genes must be expressed in at least one group (or roughly any three samples
# across the entire experiment) to be kept for downstream analysis.

library("DESeq")

keep.exprs <- rowSums(cpm>1)>=3 # is filtering is ok
DGEfiltered <- DGEobject[keep.exprs,, keep.lib.sizes = FALSE]

## Normalisation generally required to ensure that expression distributions
## of each sample are similar across the entire experiment.
## Any plot showing the per sample expression distributions, such as a density or boxplot,
## is useful in determining whether any samples are dissimilar to others.
## Normalisation by the method of trimmed mean of M-values (TMM)
## is performed using the calcNormFactors function in edgeR.
# cpmF = filtered subset # log.cpmF = log cpm filtered subset
# Copy unnormalized, filtered data to show in boxplot
DGEf.notnorm <- DGEfiltered # normaliseerde
DGEf.norm <- calcNormFactors(DGEfiltered, method = "TMM")
DGEf.norm$samples$norm.factors

# --> scaling factors relatively close to 1 --> mild TMM normalisation

## Boxplots of log-CPM values showing expression distributions
## for unnormalised data (A) and normalised data (B) for each sample
par(mfrow = c(1,1))
par(mar = c(7,4,2,1))
# Boxplot unnormalised data
DGEf.notnorm$samples$norm.factors
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
log.cpmF.notnorm <- cpm(DGEf.notnorm, log = TRUE)
boxplot(log.cpmF.notnorm, las = 2, col = sample.col, cex = 0.9,
        main = "A. Unnormalised data", ylab = "Log-cpm")
# Boxplot normalised data
DGEf.norm$samples$norm.factors

log.cpmF.norm <- cpm(DGEf.norm, log = TRUE)
boxplot(log.cpmF.norm, las = 2, col = sample.col, cex = 0.9,
        main = "B. Normalised data", ylab = "Log-cpm")
# --> only when scaling factors are not close to 1 you will see a big difference

Hello nemo, thanks for your answer, looks good. Could you add also some info where the script comes from? Maybe a reference or at least a bit of reasoning? I also believe that you should include some libraries.

I wrote it for the most part my self for a school project. The rest of the code is from the lesson I had. It was for an high throughput analysis on dataset for a disease. To determine the over expressed genes.

Cool, if the materials of the class are available it would be good to link it, but that sounds already quite good to me. I tried to [edit] your post (& add there the info), check if it still makes sense and do not hesitate to make a further edit to fix it or to add more details/links.


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

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

Then run:
Seurat::BuildClusterTree()
Seurat::FindAllMarkersNode()

Assessing the cluster markers for each node will hopefully give you a good idea on which clusters should be combined. Then you can "combine" the clusters and re-label the cells using something like:
library(plyr)
cell.labels <- seuratobject@ident
cell.labels <- mapvalues(cell.labels,
  from=0:16, # cluster numbers
  to=c('A', 'B', 'C', 'C', 'D', 'E', 'E', ... )) # etc

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

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

That is a very general recommendation. Depending on your experiment, you can get a very different number of clusters with the same number of cells at the same resolution.
You can actually use a vector of different resolutions and see which one performs best:
pbmc_small <- FindClusters(
  object = pbmc_small,
  reduction.type = "pca",
  resolution = c(0.4, 0.8, 1.2),
  dims.use = 1:10,
  save.SNN = TRUE
)

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

Have a look into clustree to assess the different clusters by clustering them and see different levels.
Example:
library(clustree)
data("iris_clusts")

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

Using the stability index to assess clusters
The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions. The stability index is automatically calculated when a clustering tree is built. Note that each level of clustering corresponds to a different resolution.
clustree(iris_clusts, prefix = "K", node_colour = "sc3_stability")


text:
Finding a single open reading frame with ribosomal binding site, using Biopython
I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a single large open reading frame contained within the DNA sequence using Biopython.
I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.
Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;
def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

OUTPUT
 TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427

However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?
My specification says to find the single open reading frame however since none of these begin the the M start codon i'm a bit lost how they can be so, I have no information on where my sequence is from (something I have to find out) so using the standard translation table 1.

From what I can tell, the program output does not indicate whether the open reading frames (ORFs) are preceded by ribosome binding sites (RBSs). I'm not sure what parameters are used when establishing ORFs. To verify your RBS requirement, index the input DNA (using given ORF locations from biopython output, pay attention to strands, which are designated by 1/-1) and search upstream for a sequence similar to your consensus. You'd have to pre-define your criteria for similarity to consensus in order to find it programmatically (or you can do it visually to better avoid false negatives). Perhaps try a search that permits ambiguous/non-matching third nucleotides in codons (like a regular expression such as "AA.GA.GT." where "." can be any nucleotide). Biologically, third nucleotide changes are least probable to drastically alter the translated residue (more likely to encode amino acid with similar chemical properties).

You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.
There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.

Thanks for the heads up, yeah I think I can safely assume that yeah the start codon will be an M. Apologies for the follow up question, is there anywhere you can point me in the right direction to in terms of achieving this? I'm struggling to find anything related.

@maltedg sorry, I just saw your comment, but I don't know what you mean. This is pointing you in the right direction. As I said, start by finding AAGGAGGTG, and then look for ORFs that start after it. If you can assume that they'll always start with an AUG, then look for AUGs close to the AAGGAGGTG and extend the reading frame from there.


text:
Downloading URL using AWK for fixed fields
I have a text file as follows
a   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz
b   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz
c   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz

The file is tab delimited and has 8 columns with 8th column being the url link  I need to download the files one by one using wget and changing the name of the downloaded file with -O option in wget. The argument I want to give to the -O option is combination col1 and col2 separated by _ so that for first line I get the downloaded file as a_k.fastq.gz and so on. For example, the first line of the file above should give:
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

Following the post from here I realized looping through the file is slow and sometimes difficult. How can I download these files using wget with awk or cut as the tool for fixed fields

Can you post up an example wget command?

wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

how can I put it into code in the comment section??

edit the original post

@Llopis we've decided that basic things that bioinformaticians need to do in their daily work are on topic even when the issue could be expressed as a general text parsing question.

Ok, I misjudged and thought it was about downloading data.

awk '{print "wget " $9 " -O " $1 "_" $2 ".fastq.gz"}' tmp

will give you
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz -O b_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz -O c_k.fastq.gz

If you add a | bash to the command it will then execute all the wget commands.

This works great. Can you tell what |bash is doing ??

@user3138373 the output of the awk command is a set of shell commands, essentially a mini script. By piping it to bash with | bash you're telling it to execute each input line as a command. So it will run the wget calls and download the sequences.

A variation of Bioathlete’s answer that doesn't suffer from escaping issues:
awk 'BEGIN{ FS="\t" } { printf("%s\t%s_%s.fastq.gz\n", $8, $1, $2) }' input |
while IFS=$'\t' read -r url filename; do
    wget -O "$filename" -- "$url"
done

Or entirely with shell parsing:
while IFS=$'\t' read -r c1 c2 c3 c4 c5 c6 c7 url; do
    wget -O "${c1}_${c2}.fastq.gz" -- "$url"
done < input

Or entirely without shell parsing but with GNU xargs:
awk 'BEGIN{ FS="\t" } { printf("%s_%s.fastq.gz\n%s\n", $1, $2, $8) }' input |
xargs -r -d '\n' -n 2 -- wget -O

… you could even run multiple wget instances in parallel from GNU xargs with the -P option, though you should set the -q option on wget to disable progress bars overwriting each other on the same terminal. Then again it would be better to use a real command-line download manager like Aria or Axel for that task.


text:
bed file to gtf conversion using awk
I have a bed file like this:
chr1    12227   12612   LOC102725121    0   +
chr1    12721   13220   LOC102725121    0   +
chr1    12227   12612   DDX11L1 0   +
chr1    12721   13220   DDX11L1 0   +
chr1    14829   14969   WASH7P  0   -
chr1    15038   15795   WASH7P  0   -
chr1    15947   16606   WASH7P  0   -
chr1    16765   16857   WASH7P  0   -
chr1    17055   17232   WASH7P  0   -
chr1    17368   17605   WASH7P  0   -

I am trying to convert into gtf file for downstream analysis. Gtf file should look like this:
chr1   myIntrons   feature 3073253   3074322   .   +   .     gene_id "intron1"; transcript_id "intron1"
chr1   myIntrons   feature 3075858   3079636   .   +   .     gene_id "intron2"; transcript_id "intron2"

For the gene_id part in the 9th column in the double quotes I want my 4th column and in transcript_id I want combination of column1,column2,and column3 so that my output looks like this:
Showing for the first line
chr1    myIntrons   intron  12228   12612   .   +   .   gene_id "LOC102725121"; transcript_id "chr1:12227-12612"

I am using this code but under double quotes it is not expanding $1 $2 etc
awk -v OFS="\t" '{print $1,"myIntrons","intron", $2+1, $3, ".", $6, ".", "gene_id $4; transcript_id $1:$2-$3"}' file.bed| head

Instead of the "gene_id ..." bit, you want sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"", $4, $1, $2, $3).


text:
how do I predict bacterial small non-coding RNA for a specific mRNA?
I am working on Vibrio parahaemolyticus. I have a gene that it is possible might be regulated by small non-coding RNAs. How can I predict possible sRNAs that target the transcript of this gene. Most software offers the opposite of what I want, i.e., they predict targets for small RNAs!
I'm trying to predict bacterial small non-coding RNA in general (Support Vector Machine) but I got more than 8 thousand sequences within intergenic regions alone. Is there a code for aligning these sequences to my transcript?

Answer from @llrs, converted from comment:
You could try to download RNAseq data from this bacteria and see if some of the RNA align to your sequence. You could try to blast your gene and look for the sequences that have low identity similarity (or classify them by type and then look what is each hit).

Good idea but from my experience, normal alignment is not optimized for sRNA. Tools like targetRNA2 are optimized for sRNA alignment such as imperfect pairing and energy ..etc.


text:
How do programs like STRUCTURE determine populations from alignments?
If I were to run BAM/SAM files through a program like STRUCTURE, or ADMIXTURE described  here how would they determine which individuals belong to what groups?
What part of the files would these programs use to determine which individuals belonged to different populations? I believe they do this via determining allele frequencies for each individual and then putting them into a population based on allele frequency, is the only criteria that they use comparison of allele frequencies between individuals?

structure works with SNP/genotype information, so won't work directly off BAM/SAM files. You would need to convert to VCF (e.g. using mpileup/vcftools), then process using something like PLINK to convert it into a format that STRUCTURE recognises.
Here's the description I have about structure from my PhD thesis. It might be a bit dated now, but I expect most of it is still applicable:
A clustering program called structure was created by Falush et al. (2003) that is used to determine population structure using a set of unlinked genotypes. The expected use case for structure is in the analysis of population structure, but it can also be useful for determining the utility of a given set of markers for categorising (or quantifying) a given trait.
The main benefit of using structure over other methods is that it is able to derive estimated ancestry coefficients from the data, rather than a discrete yes/no for each individual. This allows for the use of cutoff values for the assignment of individuals (useful for tweaking false positive and false negative rates in diagnostic tests), and a quantifiable value where a classification is not useful or does not make sense (as in the case of a continuous trait with high heritability such as height).
The program begins by assigning each individual to a given population. This can be either given in the input data as a ‘popinfo’ flag, or determined randomly by structure. After this, the program constructs a probability distribution for the data, considering the population assignments given, and uses that to estimate population allele frequencies. The program then constructs another probability distribution for these data, and uses that to estimate the population that each individual is likely to have originated from. This process is repeated many times (such a process is known as a random walk), and the underlying theory suggests that the progression of population assignments through the random walk process can be used as an approximation of the actual populations of origin for each individual.
The documentation included with structure describes the general algorithm for each iteration (of maybe hundreds of thousands of iterations) as follows:

Sample (population allele frequencies for iteration m) from the probability of (allele frequency) given (genotypes of sampled individuals) and (population of origin of individuals for iteration m − 1).
Sample (population of origin of individuals for iteration m) from the
probability of (population of origin) given (genotypes of sampled individuals) and (allele frequencies for iteration m).

In more simple terms, the steps are similar to the following statements:

these populations probably have this allele frequency distribution
therefore...
these individuals probably came from this population

The structure program is able to accommodate admixture into its model
by considering the population of origin as a probability, rather than a single
discrete value.
The file that is output from structure is a text file with individual identifiers and predicted ancestry coefficients (Q values). While the program does provide means to display these data graphically, the visual representation of data (including sorting) is better controlled using an external program. The most commonly used program is Noah Rosenberg’s distruct (Rosenberg, 2004), a program that allows colour customisation and group re-ordering. For the purpose of this thesis, a custom R script, snpchip2structure.R, has been designed which also allows selection of different sorting methods, a scatter plot for two populations, error bars, and a few additional features (see figure 1.7). More details about this script can be found in the programs directory of the supplementary CD for this thesis.


text:
Find number of possibilities
We have a column which has following format : (1, 1, 1, 0, 1).
Here is the definition of what means two columns are compatible.

Using the notation Oi to denote the collection of rows possessing a 1
  in their i-th element, we conclude that one of three possibilities
  must be true: Oi ⊆ Oj , Oj ⊆ Oi , or Oi and Oj are disjoint (the first
  two cases include the possibility that Oi = Oj). If columns i and j
  satisfy this condition, then we call them compatible

The question is : How many possible columns of length 5 are compatible with the column (1, 1, 1, 0, 1)?
You can read more here where problem is described in page 14-16.
Thanks in advance !

I'm voting to close this question as off-topic because it is not a question about bioinformatics.

@AlexReynolds, I attached a link and you will see that is about bioinformatics.

In the document there are also algorithms for finding genes tree.

That the topic it can be applied is in bioinformatics doesn't make your question about bioinformatics. You can edit explaining what is 1 and 0 and which is your bioinformatic question. But it seems a nice problem

Let $n$ be the number of rows (in your case $n=5$) and $|O_j|$ the cardinality of $O_i$.
We can calculate the number of compatible columns as sum of the number of $O_i \subseteq O_j$, $O_j \subseteq O_i$ and $O_i, O_j$ disjoint.
number of $O_i \subseteq O_j$:
Here we need to find the number of supersets of our $O_i$. That means we need to find the number of sets which contain the indices in $O_i$ and for each index not in $O_i$, we have two choices: It is contained, or it is not contained in $O_j$. Therefore the number of supersets is $2^{n-|O_i|}$, which is 2 to the power of elements not in our set.
number of $O_i, O_j$ disjoint:
Here all elements of $O_i$ are fixed not to be in $O_j$ and for each index not in $O_i$, it can again be in $O_j$ or not. That means that the number of possible sets is again $2^{n-|O_i|}$.
number of $O_j \subseteq O_i$:
Here we need to find the number of subsets of $O_i$ and this is the number of elements in its power set, which is $2^{|O_i|}$.
So finally, this means that the number of compatible columns can be calculated as: $2^{n-|O_i|} + 2^{n-|O_i|} + 2^{|O_i|} - 2$.
Why -2? We have counted 2 sets twice:

The empty set, which is a subset of $O_i$ and also counted in the number of disjoint sets
and $O_i$, which we counted in $O_i \subseteq O_j$ and in $O_j \subseteq O_i$

This means for your examle, we have $16 + 2 + 2 - 2 = 18$ compatible sets.
Note: The all zero vector (which corresponds to $O_i=\{\}$) is a special case for which the formula above does not work, as all disjunct sets and all supersets are counted twice. For this special case, the number of compatible columns is simply $2^n$.

Thanks for your answer. How do you calculate |Oi| ?

$|O_i|$ is the number of elements in $O_i$, so for your column $|O_i|$ is 4 because $O_i = {1,2,3,5}$

In my example why is not 5 ? Cardinal({1,1,1,0,1})=5

The cardinality of your column is 5. But in your definition: "Using the notation Oi to denote the collection of rows possessing a 1 in their i-th element,", $O_i$ is the set of rows with a 1. In your column there are 4 rows with a one, so the cardinality of $O_i$ = 4

Thanks a lot for your explanation.Now I understand. It means that for (0,0,0,0,0) the answer is 2^5+2^5+2^0-2 =63 ?

The number of compatible columns can not be larger than the total number of sets with 5 elements ($2^5$) :). (0,0,0,0,0) Is indeed a special case for which my rule above breaks, as all disjoint sets and all supersets are counted twice. So for This special case you have to hard code the result to $2^n$. For all other cases, the formula should work.


text:
What does it mean when a gene and its transcript have opposite orientations in a GFF3 file?
I was working with a given GFF3 file, and I observed that some transcripts have orientation opposite to their transcripts. Here is a snippet:
chr1    .       gene    2189548 2194772 .       +       .       ID=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;Name=FAAP20;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20
chr1    .       transcript      2189548 2194772 .       -       .       ID=f47feeac-74ce-1543-53b6-afd9bbb8ab19;Name=rna-NM_182533.4;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20;Parent=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;transcript_biotype=protein_coding;transcript_name=rna-NM_182533.4;protein_id=NP_872339.3

How do I interpret this? I checked UCSC Browser and in each case it seems that the browser agrees with the transcript in each case, which seems like it should be canonical.
Moreover, I observed in this file that gene orientation is always +. Whereas a highly curated GFF3 of the human genome represents both directions:
# my file
$ awk '$3 == "gene" {print}' filtered_potential_genes.gff | cut -f7 | sort | uniq -c
   1600 +
# highly curated annotations
$ awk '$3 == "gene" {print}' /MANE.GRCh38.v0.95.select_refseq_genomic.gff | cut -f7 | sort|  uniq -c
   9393 +
   9191 -

According to the person who provided me the file, the existence of bidirectional transcription was indeed the motivation for this. However, it's obviously confusing to have it still be a definitive + rather than . or something to signify unstrandedness.
It seems likely that this is due to an artifact in my file. Are there other interpretations? Is the standard simply that gene orientation is meaningless and transcript orientation is important?

Looks like a poor annotation to me. I could understand a few scattered ncRNAs being arranged in the opposite direction from the gene, but all genes being oriented identically doesn't make sense.


text:
R Biostrings pairwiseAlignment to BAM
The R package Biostrings has a function to create a pairwiseAlignment from pattern and subject sequences.
So far I can save the result into a text file using writePairwiseAlignments. I would like to save the result into a SAM/BAM file, but couldn't find anything from the BioStrings package or any R package.
What are the steps to convert the result from pairwiseAlignement into a SAM/BAM file?
Example:
library(Biostrings)

patterns = DNAStringSet(c(
  'GCTCCCTGGCGGTAAGTTGATCAAAGGAACGCAAAGTTTTCAAG',
  'ACTTCACAGCTCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG'
))

subject = DNAString('ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG')

pwa = pairwiseAlignment(patterns, subject)

writePairwiseAlignments(pwa) # produces below
#=======================================
#
# Aligned_sequences: 2
# 1: P1
# 2: S1
# Matrix: NA
# Gap_penalty: 14.0
# Extend_penalty: 4.0
#
# Length: 54
# Identity:      44/54 (81.5%)
# Similarity:    NA/54 (NA%)
# Gaps:          10/54 (18.5%)
# Score: 27.19725
#
#
#=======================================
# P1 ---------GCTCCCTGGCGGTAAGTTGATCAAAGGAA-CGCAAAGTTTTCAAG
#             ||||||||||||||||||||||||||||| |||||||||||||||
# S1 ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG

Essentially, I would like to convert the result from writePairwiseAlignments to a sam/bam format so I can view it in a GUI application. If I understand correctly, subject (referred from the Biostrings package) should be the reference, the patterns are the reads.

This really sounds like an XY problem. Why don't you just use a GUI alignment viewer? The sam format is really not designed for this.


text:
Uniprot: Determining protein and organism given a DNA sequence
I would like to find the identity of a protein, together with the organism name and kingdom.
The given DNA sequence is below:
ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTATGATATTCTTC
TTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATTCCAGATGAGTGTACCGGAGT
TAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTTAGTATGTGGTCTACTTTTGATTACGGGTTC
ATCTTGGATGGAGGAAGGTTCAGGAACGGCCTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCAT
AGCGGACTTGCTGTAGATACGTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCA
TCAACCTTATATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC
CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGCTATCACGATG
CTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGGCGGTGATCCTATGTTGTATC
AGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATATCATCATACTTCCAGTATTCGGTATAGTCAG

I've been asked to translate this DNA sequence into the 6 possible frames which are:
5'3' Frame 1
TVVR-SSVVQRPNNITCNAYDILLHHARHDGRTRQLTSAIPDECTGVSIPKD--HRYMIF
SMWSTFDYGFILDGGRFRNGLNRLSTTSAHCKS-RTCCRYVHYRIAHGRCKLPYRKHQPY
MYNRLCPPFTHGDAAVITLSLIHYNHCSVTHRSCACASRCYHDATH--KLEYQLL-QFGR
R-SYVVSALILGVWASRSLYHHTSSIRYSQ
5'3' Frame 2
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL
VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV
5'3' Frame 3
CCSVIISCTTS-QHHMQCL-YSSSSCQARWQD-ATY-CHSR-VYRS-YSQRLITSVYDF-
YVVYF-LRVHLGWRKVQERPEPSIHH-RSLQVIADLL-IRSLSHCTWPVQAPLQEASTLY
VQSPMPAVHSWRCCSHHFIPDPLQSLQRYS-ELCLC-QVLSRCYSLIEVGVPASMTVRQA
VILCCISTYSGCLGIQKSISSYFQYSV-S
3'5' Frame 1
LTIPNTGSMMI-TSGCPNTQNKC-YNIGSPPAELS-KLVLQLLSVSSIVIAPASTGTTPM
SNAAVIVMDQG-SDDCSIAMSERRA-AIVHIRLMLPVRELAPAMCNAIMNVSTASPL-LA
VSASGG-TVQAVPEPSSIQDEPVIKSRPHTKNHIPMLLIFGNTNSGTLIWNGTSKLPSPA
IVPGMMKKNIISIACDVVRTLYN--SPNNS
3'5' Frame 2
-LYRILEV--YRLLDAQTPRISADTT-DHRLPNCHRSWYSNFYQ-VAS--HLLAQAQLL-
VTLQ-L-WIRDKVMTAASP-VNGGHRRLYI-G-CFL-GSLHRPCAMR--TYLQQVRYDLQ
-ALVVDRRFRPFLNLPPSKMNP-SKVDHILKIIYRCY-SLGILTPVHSSGMALVSCLVLP
SCLA--RRIS-ALHVMLLGRCTTDDHRTT
3'5' Frame 3
DYTEYWKYDDIDFWMPKHPE-VLIQHRITACRTVIEAGTPTSISE-HRDSTC-HRHNSYE
-RCSDCNGSGIK--LQHRHE-TAGIGDCTYKVDASCKGACTGHVQCDNERIYSKSAMTCS
ER-WWIDGSGRS-TFLHPR-TRNQK-TTY-KSYTDVINLWEY-LRYTHLEWH--VA-SCH
RAWHDEEEYHKHCM-CC-DVVQLMITEQQ

I ran jobs on Uniprot for all these sequences. Here are the results:
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E0FZ
www.uniprot.org/blast/uniprot/B2017100383C3DD8CE55183C76102DC5D3A26728BEBF5BAE
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E5FP
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082F08D
nothing
www.uniprot.org/blast/uniprot/B20171003F725F458AC8690F874DD868E4ED79B88CE32C4B

But how can I decide the identity of the protein and the name of the organism?
In the first one, the organism is Pacific ridley turtle  but in the second one it is Diplonema sp. ATCC 50224. So I can only decide at a kingdom level, but not to the organism, right?

Hi snnlaynnkrdsm, I've edited your question to make it sound less "homeworkey", as I suspect that would put people off answering it.

To be clear, homework questions are fine on this site. But you should be clear when something is a homework assignment. Accordingly, we often won't provide answers to homework questions outright, but lead you in the right direction.

Well, let's see. The first step is to simply run a nucleotide blast (blastn) against the nr (non-redundant) database. I did this using the NCBI blast tool. While, as expected, I got many hits, there was only one hit which matched the entire length of the query sequence perfectly:
Diplonema sp. ATCC 50224 cytochrome c oxidase subunit 1 (cox1) mRNA, complete cds; mitochondrial
Sequence ID: JF698657.1
Length: 1612Number of Matches: 1
Range 1: 141 to 770
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
1164 bits(630)  0.0 630/630(100%)   0/630(0%)   Plus/Plus
Query  1    ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  60
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  141  ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  200

Query  61   GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  120
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  201  GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  260

Query  121  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  180
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  261  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  320

Query  181  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  240
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  321  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  380

Query  241  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  300
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  381  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  440

Query  301  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  360
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  441  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  500

Query  361  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  420
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  501  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  560

Query  421  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  480
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  561  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  620

Query  481  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  540
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  621  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  680

Query  541  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  600
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  681  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  740

Query  601  CATCATACTTCCAGTATTCGGTATAGTCAG  630
            ||||||||||||||||||||||||||||||
Sbjct  741  CATCATACTTCCAGTATTCGGTATAGTCAG  770

So, we now know the sequence is from the gene for mitochondrial Cytochrome C Oxidase subunit 1 from Diplonema sp. ATCC 50224. That the sequence is mitochondrial, is very important since the mitochondrion of this species has a slightly different genetic code. If you go to the NCBI's page for this sequence (JF698657.1), you will see this line:
/transl_table=4

That tells you it is using translation table 4:
4. The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code (transl_table=4)

    AAs  = FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
  Starts = --MM------**-------M------------MMMM---------------M------------
  Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG

See how all of your translated reading frames have in-frame stop codons? That's because you translated them using the standard genetic code instead of the correct mitochondrial one. As you can see above, in this code, the codon TGA is not a STOP, as in the standard code, and instead codes for tyrosine (W).
So, if you look at the result of the protein blast you ran on UniProt, you will see:
Query1   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL   59
         LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI FL
F5CI61   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL   107

Query60  VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   118
         VCGLLLITGSSWMEEGSGTA TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
F5CI61   VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   167

Query119 CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   177
         CTIAYARRSLMAMLQSSLYP SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
F5CI61   CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   227

Query178 GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  206
         GDPMLYQHLFWVFGHPEVYIIILPVFGIV
F5CI61   GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  256

The only locations where there isn't a 100% identical match is where a W is aligned to a STOP (- in the Uniprot output). So, your actual protein sequence is the 2nd frame of the 5' to 3' direction which, when translated using the correct code, is:
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL
VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV

Well if it is not a homework assignment, I would suggest to use blastx. Then you don't need to translate your DNA into protein first, but that will be covered by blastx.
If it is a homework assignment though, the teacher probably wants you to select the best blast hit. So which hit gives the highest score, and most matching (identical) amino acids between the alignment.


text:
How can I extract gene names for a metabolic pathway from KEGG?
Note: this question has also been asked on Biostars
I need the get the list of gene names involved in glycolysis (to put an example). Not manually, I need to do this in a script. Ideally with Python. I think KEGG is the proper database to do this, but any other (free) suggestion would do.
How can I do it?

To your knowledge there are other metabolic pathways, including metaCyc, reactome,...

I could see that this question is reposted on Biostars. Here is the solution, https://www.biostars.org/p/274321/#274326

How to do vice verse, if i want to find all pathways associated with a list of genes using R or python? Please advise

using the REST API of togows: one can fetch a json structure of the genes associated to a kegg pathway. eg: http://togows.org/entry/kegg-pathway/hsa00010/genes.json
[
  {
    "3101": "HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]",
    "3098": "HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]",
    "3099": "HK2; hexokinase 2 [KO:K00844] [EC:2.7.1.1]",
    "80201": "HKDC1; hexokinase domain containing 1 [KO:K00844] [EC:2.7.1.1]",
    "2645": "GCK; glucokinase [KO:K12407] [EC:2.7.1.2]",
    "2821": "GPI; glucose-6-phosphate isomerase [KO:K01810] [EC:5.3.1.9]",
    "5213": "PFKM; phosphofructokinase, muscle [KO:K00850] [EC:2.7.1.11]",
    "5214": "PFKP; phosphofructokinase, platelet [KO:K00850] [EC:2.7.1.11]",
    "5211": "PFKL; phosphofructokinase, liver type [KO:K00850] [EC:2.7.1.11]",
    "2203": "FBP1; fructose-bisphosphatase 1 [KO:K03841] [EC:3.1.3.11]",
    "8789": "FBP2; fructose-bisphosphatase 2 [KO:K03841] [EC:3.1.3.11]",
    "230": "ALDOC; aldolase, fructose-bisphosphate C [KO:K01623] [EC:4.1.2.13]",
    "226": "ALDOA; aldolase, fructose-bisphosphate A [KO:K01623] [EC:4.1.2.13]",
    "229": "ALDOB; aldolase, fructose-bisphosphate B [KO:K01623] [EC:4.1.2.13]",
    "7167": "TPI1; triosephosphate isomerase 1 [KO:K01803] [EC:5.3.1.1]",
    "2597": "GAPDH; glyceraldehyde-3-phosphate dehydrogenase [KO:K00134] [EC:1.2.1.12]",
    "26330": "GAPDHS; glyceraldehyde-3-phosphate dehydrogenase, spermatogenic [KO:K10705] [EC:1.2.1.12]",
    "5232": "PGK2; phosphoglycerate kinase 2 [KO:K00927] [EC:2.7.2.3]",
    "5230": "PGK1; phosphoglycerate kinase 1 [KO:K00927] [EC:2.7.2.3]",
    "5223": "PGAM1; phosphoglycerate mutase 1 [KO:K01834] [EC:5.4.2.11]",
    "5224": "PGAM2; phosphoglycerate mutase 2 [KO:K01834] [EC:5.4.2.11]",
    "441531": "PGAM4; phosphoglycerate mutase family member 4 [KO:K01834] [EC:5.4.2.11]",
    "2027": "ENO3; enolase 3 [KO:K01689] [EC:4.2.1.11]",
    "2026": "ENO2; enolase 2 [KO:K01689] [EC:4.2.1.11]",
    "2023": "ENO1; enolase 1 [KO:K01689] [EC:4.2.1.11]",
    "5315": "PKM; pyruvate kinase, muscle [KO:K00873] [EC:2.7.1.40]",
    "5313": "PKLR; pyruvate kinase, liver and RBC [KO:K12406] [EC:2.7.1.40]",
    "5161": "PDHA2; pyruvate dehydrogenase alpha 2 [KO:K00161] [EC:1.2.4.1]",
    "5160": "PDHA1; pyruvate dehydrogenase alpha 1 [KO:K00161] [EC:1.2.4.1]",
    "5162": "PDHB; pyruvate dehydrogenase (lipoamide) beta [KO:K00162] [EC:1.2.4.1]",
    "1737": "DLAT; dihydrolipoamide S-acetyltransferase [KO:K00627] [EC:2.3.1.12]",
    "1738": "DLD; dihydrolipoamide dehydrogenase [KO:K00382] [EC:1.8.1.4]",
    "160287": "LDHAL6A; lactate dehydrogenase A like 6A [KO:K00016] [EC:1.1.1.27]",
    "92483": "LDHAL6B; lactate dehydrogenase A like 6B [KO:K00016] [EC:1.1.1.27]",
    "3939": "LDHA; lactate dehydrogenase A [KO:K00016] [EC:1.1.1.27]",
    "3945": "LDHB; lactate dehydrogenase B [KO:K00016] [EC:1.1.1.27]",
    "3948": "LDHC; lactate dehydrogenase C [KO:K00016] [EC:1.1.1.27]",
    "124": "ADH1A; alcohol dehydrogenase 1A (class I), alpha polypeptide [KO:K13951] [EC:1.1.1.1]",
    "125": "ADH1B; alcohol dehydrogenase 1B (class I), beta polypeptide [KO:K13951] [EC:1.1.1.1]",
    "126": "ADH1C; alcohol dehydrogenase 1C (class I), gamma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "131": "ADH7; alcohol dehydrogenase 7 (class IV), mu or sigma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "127": "ADH4; alcohol dehydrogenase 4 (class II), pi polypeptide [KO:K13980] [EC:1.1.1.1]",
    "128": "ADH5; alcohol dehydrogenase 5 (class III), chi polypeptide [KO:K00121] [EC:1.1.1.1 1.1.1.284]",
    "130": "ADH6; alcohol dehydrogenase 6 (class V) [KO:K13952] [EC:1.1.1.1]",
    "10327": "AKR1A1; aldo-keto reductase family 1 member A1 [KO:K00002] [EC:1.1.1.2]",
    "217": "ALDH2; aldehyde dehydrogenase 2 family (mitochondrial) [KO:K00128] [EC:1.2.1.3]",
    "224": "ALDH3A2; aldehyde dehydrogenase 3 family member A2 [KO:K00128] [EC:1.2.1.3]",
    "219": "ALDH1B1; aldehyde dehydrogenase 1 family member B1 [KO:K00128] [EC:1.2.1.3]",
    "501": "ALDH7A1; aldehyde dehydrogenase 7 family member A1 [KO:K14085] [EC:1.2.1.3 1.2.1.8 1.2.1.31]",
    "223": "ALDH9A1; aldehyde dehydrogenase 9 family member A1 [KO:K00149] [EC:1.2.1.3 1.2.1.47]",
    "221": "ALDH3B1; aldehyde dehydrogenase 3 family member B1 [KO:K00129] [EC:1.2.1.5]",
    "222": "ALDH3B2; aldehyde dehydrogenase 3 family member B2 [KO:K00129] [EC:1.2.1.5]",
    "220": "ALDH1A3; aldehyde dehydrogenase 1 family member A3 [KO:K00129] [EC:1.2.1.5]",
    "218": "ALDH3A1; aldehyde dehydrogenase 3 family member A1 [KO:K00129] [EC:1.2.1.5]",
    "84532": "ACSS1; acyl-CoA synthetase short chain family member 1 [KO:K01895] [EC:6.2.1.1]",
    "55902": "ACSS2; acyl-CoA synthetase short chain family member 2 [KO:K01895] [EC:6.2.1.1]",
    "130589": "GALM; galactose mutarotase [KO:K01785] [EC:5.1.3.3]",
    "5236": "PGM1; phosphoglucomutase 1 [KO:K01835] [EC:5.4.2.2]",
    "55276": "PGM2; phosphoglucomutase 2 [KO:K15779] [EC:5.4.2.7 5.4.2.2]",
    "2538": "G6PC; glucose-6-phosphatase catalytic subunit [KO:K01084] [EC:3.1.3.9]",
    "57818": "G6PC2; glucose-6-phosphatase catalytic subunit 2 [KO:K01084] [EC:3.1.3.9]",
    "92579": "G6PC3; glucose-6-phosphatase catalytic subunit 3 [KO:K01084] [EC:3.1.3.9]",
    "83440": "ADPGK; ADP dependent glucokinase [KO:K08074] [EC:2.7.1.147]",
    "669": "BPGM; bisphosphoglycerate mutase [KO:K01837] [EC:5.4.2.11 5.4.2.4]",
    "9562": "MINPP1; multiple inositol-polyphosphate phosphatase 1 [KO:K03103] [EC:3.1.3.80 3.1.3.62]",
    "5105": "PCK1; phosphoenolpyruvate carboxykinase 1 [KO:K01596] [EC:4.1.1.32]",
    "5106": "PCK2; phosphoenolpyruvate carboxykinase 2, mitochondrial [KO:K01596] [EC:4.1.1.32]"
  }
]

I whipped this little script up using the KEGG API:
#!/usr/bin/env python3
import urllib.request
import re
import sys
pathway = 'hsa00010' # glycolysis
url = "http://rest.kegg.jp/get/" + pathway
with urllib.request.urlopen(url) as f:
    lines = f.read().decode('utf-8').splitlines()
    want = 0
    for line in lines:
        fields = line.split()
        ## The list of genes starts here
        if fields[0] == 'GENE':
            want = 1
            ## The line with GENE is different
            print(fields[2].rstrip(';'))
        ## We reached the next section of the file
        elif want == 1 and re.match('^\S', line):
            sys.exit();
        ## We're still in the list of genes
        if want == 1 and len(fields)>1:
            print(fields[1].rstrip(';'))

Answers from Biostars:
GeneSCF 'prepare_database' module will extract all the pathways with corresponding genes as simple table format in plain text file. Check 'Two step process' (sub-heading) from the provided link (First step will do your job).

You can use KEGG REST API for extracting gene names from the pathway. Python request library can be used to make the REST calls to KEGG REST API and query information. Please follow this tutorial if you are a newbie to the Python programming and REST web services.


text:
pvalue calculation of protein-protein network with permutation test
Note: this question has also been asked on Biostars
I have performed a network recreation analysis based on the interactions of proteins from String db.
I wanted to know if the interactions have significance by reshuffling the proteins with permutation against the background network. Since, its first time for me to perform permutation test on PPI, Can someone guide me to a tool or a tutorial to perform permutation test on newly created PPI network.
I tried looking at a youtube tutorial and read some articles like this one, but they were not very useful for my purpose.

What kind of permutation test you want to do with your network? What do you want to permutate (and in which programming language)? In R you could use boost to permute

Please [edit] your question and explain exactly what type of permutation you need. Do you need to keep the edges the same and randomize the nodes? The nodes the same and the edges randomized? Randomize both but keep the same number of nodes/edges? Something else?

There is probably not a tool so you'll have to code it yourself. You need to provide more information to help everyone help you. How did you come across this list of proteins? Having a list of differentially expressed proteins that are known to interact in STRING might be enough to prove your hypothesis. Do you have an appreciation for when permutations are appropriate? Will doing permutations help prove or disprove your hypothesis? What do you mean by reshuffling the proteins? See Terdon's comment. Are you permuting the proteins in your list and not the sample labels?

It is probably necessary to correct for the degree distribution of the network, I am not sure that a tool such as the Mantel test will quite do this (as suggested by user Michael G.). In naturally occurring networks degree distribution is a strong confounder.
Then again I am not totally sure what shape your data is in or what the question is. You mention a background network, but I am not sure what this is or how you want to use it. Possibly your network reconstruction is a subset of your background network? In which case you just define the rule for how your subsample is drawn and then you compare it to permuted networks.
However, here is some fairly crude but serviceable R code that implements a general purpose network rewiring test, using igraph. I am not using a background network as again I'm not sure what it is.
(Python also has igraph so should be similar rewire() function)
compute_measure = function(network) {
  # i don't know what you are measuring, so you have to write this function.
  # this function takes a graph object as input
  # this function returns a single numeric value (whatever you're measuring) as output
}

# "net" is a network of the graph type from igraph

measure = compute_measure(net)

num_perms = 1000

perm_measures = c()

for (i in 1:num_perms) {
  # permute the network by rewiring, rewire 10x as many times as edges
  # to ensure saturation
  perm_net = rewire(net, with = keeping_degseq(niter = gsize(net) * 10))
  perm_measures[i] = compute_measure(perm_net)
}

# one-sided monte carlo p value is proportion of rewirings (plus one) greater than your measure
p_value = (1 + length(which(perm_measures >= measure))) / num_perms

If on the other hand you are just trying to come up with a single permutation of the edges of your network to compare against, you can run rewire() only once to get that permuted network. It depends on what question you are trying to answer, and what null hypothesis you're trying to get a p value for.

Two answers:
What you have asked for is a straight Mantel's test. You could read up on the test, it is traditionally used to assess genetic vs geographic distance but could be applied.
This is easy to do in R:
mantel.test(m1, m2, nperm = 999, graph = FALSE,
            alternative = "two.sided",  ...)

You will need to import your data (recommend csv) and parse to a data.frame both for your observed (m1) and expected (m2). Its easy.
The documentation is here.
My advice is to think carefully about the contruction of the matrix: its an everything versus everything comparison.
Alternative approach
My personal choice would be to use a covariance/correlation matrix. That can be done in R too, but you need to understand what covariance is and how it relates to your problem. If you don't really get covariance, I'd leave it. Hence I'll avoid putting code down - but is part of ggplot2.
Please note, I don't supply detailed codes (I'm not your supervisor), unless its Perl, which this isn't


text:
Sequence alignment using BWT
My Problem:
Skipping some specific background, what I want to do is judging whether some soft-clipping sequences are the same, which may result by the same SV event. Colored bases in Fig.1 is an example of soft-clipping sequences.

I use BWA -MEM as the aligner. And I am trying to develop an SV caller actually, what I want to do here is refining PE(Paired-Ends) signal. If it is SR(Split-Reads), I don't need to do this, because of SR and its supplemental alignments coming from the same read. But PE signal generated from some cluster procedures need some filters. And I think K-mer is not a good way to do this which adapted by some SV callers. For me, as a bioinformatic fool, imitation and implementation may be a good learning process. It seems still a long way to go.
My solution:
I suppose this is a sequence alignment problem. Then I use FM-indexing based on BWT to do this. My inexact/approximate matching implementation is referred from BWA.
My workflow:

Extract these soft-clipping sequences
Choose the longest one as the reference
Iterate others to match with the reference, if the number of unmatched is large than 1, it means these sequences are not the same.

Question 1:
Is that sounds reasonable to you?  Any suggestions?

Question 2:
If it is OK, how should I implement the penalty scoring scheme?
Here, because these sequences should come from the same SV event, the main differences(Mismatch, InDel) should arise from sequencing errors, not small genomic variants. So I want to set a higher gap open penalty and even higher gap extend penalty, like
mismatch_penalty = 4
gapopen_penalty = 11
gapextend_penalty = 20
alignment penalty cutoff = 30

The later means that one alignment with two more gaps will be considered as spurious alignment.
How do I optimize these parameters?

You just need to check the flag for suplementary mapping. If there is suplementary mapping then use the coordinates; if they are all the same the seq is the same and you just have one event, if different reads have different suplementary mapping locations then there's multiple events. The problem is more complex than this though so just use an existing tool.


text:
Help with generating network from gene expression data
I'm struggling to figure out how to go from gene expression data to  be able to create a network out of the most expressed genes in a cell line. Ok, I have developed an application for Cytoscape where I want to compare my analysis with research on the NCI-60 cell lines. My goal is to do these steps:
1: Use one of the cell lines in the NCI-60 cell lines, in my case I will use the NCI_H23 where there is 3 samples in this dataset.
2: Use the most expressed genes in the NCI_H23 cell line to create a subnetwork using a protein interaction network database like iRefIndex, etc.
After these two steps I should finally be able to upload the network to Cytoscape and perform my analysis, but being a programmer  I'm kind of stuck. I have tried using one of the 3 samples where I mapped the Affymatrix IDs to gene symbol and then selecting the top 20% genes with the highest values using a cutoff value, but I'm not sure if this is the correct way to do it and I'm not sure how to do step nr.2 above if my list of genes is correct.
My application investigates drug combinations based on graph theory so I want to create a network that represents a cell-line so I can compare my scoring of combinations with the ones from the NCI Almanac study. I don't want to compare it with other cell-lines, but rather find the most expressed genes in that cell-line to then retrieve interaction data from iRefIndex to get edges for the network. I guess I should average the 3 samples first before I repeat what I listed above.
I would really appreciate help on this as I've been stuck on this for over a week now.
What I'm stuck at now is trying to create a network from iRefIndex (using iRefR) using the set of gene symbols that I have extracted from the 3 samples. To be honest, I'm finding it hard to find the right tools and assessing what kind of interaction data i should use from iRefIndex etc, so I guess I am looking for confirmation that what I'm trying makes sense.

You seem to make quite a lot of guesses, maybe you should seek guidance with an advisor (offline). Usually one selects the genes that have a change in the expression when comparing with and without the drug. Have you read how other people investigate drug combinations? I don't know why have you selected the iRefIndex database, but are you looking for a interaction database or a pathway database?


text:
Help with identifying disease modules
I've made an application that at this point ranks all combinations of drug pairs relevant to a biological network/graph in the order of how disruptive the outcome of deleting the targets of a given drug pair is. So, in short, my applications objective is to kill/disrupt the network, but the optimal goal would be to disrupt nodes related to disease and have minimal impact on the functional part of the cell/network. I've read a few articles about network medicine and found that a diffusion-based method for identifying disease modules seems easier to understand coming from an informatics background.
I know how to implement the algorithm for random walkers with reset algorithm (RWR) but unsure how to get/ identify known disease associated genes/nodes of the network and these are used as seeds for the algorithm.
If you guys have some advice on how to do this or even know of a better way to identify disease, functional and topological modules I would appreciate it.
Some details of my project:
Application developed to run in Cytoscape, programmed in Java. Using Jung and JGraphT for graph algorithms.
This is for exploring paired drug combinations for treating cancer. I retrieve drug-target relation from a refined list published in an article called: "A comprehensive map of molecular drug targets" from December.


text:
How to proceed after contigs.fa file obtained from VELVET assembly
I want to assemble a linear genome and I have a reference genome too.
I ran the following commands for de novo velvet assembly:
velveth auto 31,45,2 -fastq -shortPaired1 -separate 2-CA-SVA-2_S2_L001_R1_001.fastq 2-CA-SVA-2_S2_L001_R2_001.fastq

I chose auto_33 for further analysis
velvetg auto_33 -exp_cov auto

Now my output files inside auto_33 folder are:
contigs.fa
Graph2
LastGraph
Log
PreGraph
Roadmaps
Sequences
stats.txt

But I have no clue how to move ahead. Can I do anything else with velvet or do I have to switch to another program?

Blat the contigs

Would you mind elaborating ?

How to proceed depends entirely on  what you want the assembly for. Are you looking for variants? Is this the genome of a related species which has not been annotated and you want to find its genes? Are you looking for one gene in particular? For syntenic regions? We can't tell you how to proceed if you don't tell us what you need. Please [edit] your question and clarify.

The contigs.fa file is a FASTA formatted file of the final contigs generated by velvet.  I would start assessing the quality of the assembly that was generated. Here are some questions to get you started.  Really the next steps depend on how you plan to use the assembly.
What is the distribution of contig lengths? Average and N50 contig length
How does the assembly compare to the reference?

How much of the reference does the assembly cover?

When you align the original reads back to the contigs how well are they covered?

Do any of the read alignments not make sense?
Are there sections of the contigs that are not covered?


text:
Fastest way to demultiplex bam file based on field
I have a bam file with aligned reads from multiple experiments. For each experiment I have a portion of the sequence identifier, thus for example

HWI-D00310:448:CCG0KANXX:7:1101:2570:1976,

and is different from other the other of single-run flow cells in the analysis.
I want to demultiplex the input bam file into different bam files depending on the experiment. My only concern is a fast methodology to demultiplex a bam using combinations of single-run, flow-cell, and lane. What I came up with (the fastest) is something like:
# Parallelized per chromosome. Assuming only 2 plates here.
printf '%s\n' $(seq 1 22) X Y MT \
| parallel 'samtools view input.bam {} \
  | tee >(grep -f patterns_experiment1.txt > experiment1_chr{}.sam.noheader) \
        >(grep -f patterns_experiment2.txt > experiment2_chr{}.sam.noheader)' \
&& seq 1 2 \
| parallel 'samtools view -H input.bam && cat experiment{}_*sam.noheader \
  | samtools view -1 - > input_experiment{}.bam'

Is there a better/faster way to do it?
If I understood correctly read groups should be specified prior to mapping? I started from the already mapped BAM unfortunately, but looking at the original FASTQs I noticed the pattern in which some parts of the sequence identifier were unique for one experiment vs the others.


text:
Given a multiple sequence alignment, how do I output all of the non-consensus characters by location?
I have an MSA from MAFFT (as a FASTA file; there are several sequences in the MSA), and I would like to quantify the "variants" between these sequences.
What is the standard tool to do something like this?
I would prefer to know the non-consensus variants in the the entire MSA, but that may be too tricky, so I would like to figure out the non-consensus variants between those select sequences.
I'm familiar with calculating an MSA, and then visualizing it---but I would like an output file with all of the different bases.

This is surprisingly not done or requested in the community, I'm also struggling.

A very basic solution in python:
from Bio import SeqIO
seq_dict = {i : rec.seq for i,rec in enumerate(SeqIO.parse("your_msa.aln", "fasta"))}
for pos in range(len(seq_dict[0])):
    pos_bases = ([seqio[i] for seqio in [*seq_dict.values()]])
    # If position contains more than one unique base
    if len([*set(pos_bases)])>1:
        print(pos, pos_bases)        # prints the position and alternative bases there

If you need something more useful, just leave a comment here and I'll develop this further.


text:
Trying to show the gaps of each seq on bio::Graphics after converting clustalw
I want a box representing each sequence, positioned as they are in the alignment and with gaps shown as breaks in the each box.
I've been having trouble for a while with this and have been trying to get this on my own, but I can't. I'm trying to display the gaps using Bio::Graphics after retrieving my sequences from an alignment file.
I removed the exit portion. I was trying to format the code for this website and accidentally put the exit there. The input fill is a custalw file which is an alignment file with sequences on it. The sequences have "-" that represent gaps.
The command is ./aln.pl test-bioaln.aln > ll.png
My expected output would be a .png file with sequences that are aligned and show gaps as a different color.
The input file is:
CLUSTAL W (1.81) multiple sequence alignment

JD1:1:102:1601:ORFJ00027      ------------------------------atgtataaacaacaatattttatttct--c
94a:1:107:117:orf00001        ------------------------------atgtataaacaacaatattttatttct-ac
118a:1:106:158122218:orf00020 ------------------------------atgtataaacaacaatattttatttct-gc
B31:1:100:4091:ORFB0018       ------------------------------atgtataaacaacaatattttatttctggc
72a:1:105:32:orf00022         ------------------------------atgtataaacaacaatattttatttctggc
64b:1:110:473:orf00001        ------------------------------atgtataaacaacaatattttatttctggc
29805:1:108:171:orf00001      ------------------------------atgtataaacaacaatattttatttctggc
BOL26:1:111:60:orf00001       ------------------------------atgtataaacaacaatattttatttctggc
CA-11.2A:1:109:33:orf00001    ------------------------------atgtataaacaacaatattttatttctggc
WI91-23:1:112:493:orf00001    ------------------------------atgtataaacaacaatattttatttctggc
297:1:103:411:ORFB00012       ttggatagattttatacaaagaaggtaataatgtataaacaacaatattttatttctggc
N40:1:101:1716:ORFK00021      ------------------------------atgtataaacaacaatattttatttctggc
ZS7:1:113:22:orf00001         ------------------------------atgtataaacaacaatattttatttctggc
                                                            ******************************

JD1:1:102:1601:ORFJ00027      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
94a:1:107:117:orf00001        aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
118a:1:106:158122218:orf00020 aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
B31:1:100:4091:ORFB0018       aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
72a:1:105:32:orf00022         aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
64b:1:110:473:orf00001        aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
29805:1:108:171:orf00001      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
BOL26:1:111:60:orf00001       aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
CA-11.2A:1:109:33:orf00001    aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
WI91-23:1:112:493:orf00001    aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
297:1:103:411:ORFB00012       aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
N40:1:101:1716:ORFK00021      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
ZS7:1:113:22:orf00001         aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
                              ***** *********************** ******************************

JD1:1:102:1601:ORFJ00027      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
94a:1:107:117:orf00001        aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
118a:1:106:158122218:orf00020 aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
B31:1:100:4091:ORFB0018       aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
72a:1:105:32:orf00022         aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
64b:1:110:473:orf00001        aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
29805:1:108:171:orf00001      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
BOL26:1:111:60:orf00001       aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
CA-11.2A:1:109:33:orf00001    aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
WI91-23:1:112:493:orf00001    aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
297:1:103:411:ORFB00012       aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
N40:1:101:1716:ORFK00021      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
ZS7:1:113:22:orf00001         aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
                              ******************** ***************************************

JD1:1:102:1601:ORFJ00027      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
94a:1:107:117:orf00001        aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
118a:1:106:158122218:orf00020 aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
B31:1:100:4091:ORFB0018       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
72a:1:105:32:orf00022         aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
64b:1:110:473:orf00001        aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
29805:1:108:171:orf00001      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
BOL26:1:111:60:orf00001       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
CA-11.2A:1:109:33:orf00001    aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
WI91-23:1:112:493:orf00001    aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
297:1:103:411:ORFB00012       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
N40:1:101:1716:ORFK00021      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
ZS7:1:113:22:orf00001         aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
                              ********************************************************* **

JD1:1:102:1601:ORFJ00027      aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
94a:1:107:117:orf00001        aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
118a:1:106:158122218:orf00020 aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
B31:1:100:4091:ORFB0018       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
72a:1:105:32:orf00022         aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
64b:1:110:473:orf00001        aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
29805:1:108:171:orf00001      aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
BOL26:1:111:60:orf00001       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
CA-11.2A:1:109:33:orf00001    aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
WI91-23:1:112:493:orf00001    aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
297:1:103:411:ORFB00012       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
N40:1:101:1716:ORFK00021      aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
ZS7:1:113:22:orf00001         aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
                              **************** *******************************************

JD1:1:102:1601:ORFJ00027      tattat
94a:1:107:117:orf00001        tattat
118a:1:106:158122218:orf00020 tattat
B31:1:100:4091:ORFB0018       tattat
72a:1:105:32:orf00022         tattat
64b:1:110:473:orf00001        tattat
29805:1:108:171:orf00001      tattat
BOL26:1:111:60:orf00001       tattat
CA-11.2A:1:109:33:orf00001    tattat
WI91-23:1:112:493:orf00001    tattat
297:1:103:411:ORFB00012       tattat
N40:1:101:1716:ORFK00021      tattat
ZS7:1:113:22:orf00001         tattat
                              ******

Here's the code:
#!/usr/bin/perl
use Bio::AlignIO;
use Bio::Graphics::Panel;

my $line = shift @ARGV;
my $in = Bio::AlignIO->new(-file=>$line,-format=>"clustalw");
while($aln = $in->next_aln()){
    foreach $seqobj($aln->each_seq()){
        my $seq = $seqobj->seq;
        my $id = $seqobj->id;
        my $length = $seqobj->length;
        my $seqobj = Bio::SeqFeature::Generic->new(-start =>1, -end=>$length,-display_name=>$id);
        push (@seq, $seqobj);
    }
    foreach $seq(@seq){
        my @features = $seq->get_SeqFeatures;
        my %sorted_features;
        for my $f (@features) {
            my $tag = "\-";
            push @{$sorted_features{$tag}},$f;
        }

        my $panel = Bio::Graphics::Panel->new(
            -length    => $seq->length,
            -key_style => 'between',
            -width     => 800,
            -pad_left  => 10,
            -pad_right => 10,
        );
        $panel->add_track(generic => Bio::SeqFeature::Generic->new(-start=>1,
                -end=>$seq->length),
            -glyph  => 'generic',
            -bgcolor => 'blue',
            -label  => 1,
        );

        my @colors = qw(cyan orange blue purple green chartreuse magenta yellow aqua);
        my $idx    = 0;
        for my $tag (sort keys my %sorted_features) {
            my $features = $sorted_features{$tag};
            $panel->add_track($features,
                -glyph    =>  'generic',
                -bgcolor  =>  $colors[$idx++ % @colors],
                -fgcolor  => 'black',
                -font2color => 'red',
                -key      => "${tag}s",
                -bump     => +1,
                -height   => 8,
                -label    => 1,
                -description => 1,
            );
        }
        print $panel->png;

    }
}

Unfortunately, $seq->get_SeqFeatures never returns anything. I'm not sure why it doesn't since there should be features on it coming from the Bio::SimpleAlign object. Maybe I should set some new parameters?