text:
How can I improve ESS values for ucldMean in beast?
Using starBeast I ran a 50 locus tree, with 34 individuals, using evolutionary models indicated by AICc in jModel test. I ran this tree with a MCMC chain length of 500 million, storing every 10,000 generations. I used a relaxed clock log normal.
I get low ESS values for posterior, prior, birthRate.t:Species, YuleModel.t:Species, popMean, and TreeHeight.Species.
ucldMean.c:Locus_name is low for each loci.
The TreeHeight statisitcs also have low ESS values.
Also rate.c:Locus_name.mean for each loci is low as well.
Is there a way to change my parameters in Beauti to help increase the ESS values of these statistics? Is there a parameter that would likely cause the ESS values of these particular statistics to be low? Would increasing the MCMC chain length increase the ESS values, given that I've already run half a billion?

It is not my field, but I don't understand many abbreviations you use. You might get (more) answers if you explain with simpler words what is the problem. As I understand you want to increase the value of one value by optimizing some parameters of a program, is this helpful? Is there something in the mailing list of the tool helpful? Did you try or read something to solve this?

The official BEAST website suggests a number of solutions to improve the ESS of a parameter
It is not uncommon to have several billions of generations so I don't imagine increasing your number of generations would hurt.


text:
Hisat2 : which option should mention for strand specific library read
I was trying hisat2 I get confused by the strand options...
The question has been already asked on github here but didn't get any satisfactory answer.
I'm surprise how much this strand information could be confusing in different tools.
Looking at the hisat2 manual I get confuse by the two different options:
--rna-strandness

Specify strand-specific information: the default is unstranded.
  For single-end reads, use F or R. 'F' means a read corresponds to a transcript. 'R' means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF.
  With this option being used, every read alignment will have an XS attribute tag: '+' means a read belongs to a transcript on '+' strand of genome. '-' means a read belongs to a transcript on '-' strand of genome.

and --fr/--rf/--ff

The upstream/downstream mate orientations for a valid paired-end alignment against the forward reference strand. E.g., if --fr is specified and there is a candidate paired-end alignment where mate 1 appears upstream of the reverse complement of mate 2 and the fragment length constraints (-I and -X) are met, that alignment is valid. Also, if mate 2 appears upstream of the reverse complement of mate 1 and all other constraints are met, that too is valid. --rf likewise requires that an upstream mate1 be reverse-complemented and a downstream mate2 be forward-oriented. --ff requires both an upstream mate 1 and a downstream mate 2 to be forward-oriented. Default: --fr (appropriate for Illumina's Paired-end Sequencing Assay).

Should both need to be specified ? Are they redundant ?
I didn't find so many person struggling with that, so I guess there is something obvious that I don't get. I hope someone may enlighten me about it !

--fr/--rf/--ff should rarely be set, since they refer to the relative orientation of reads and basically everything these days is --fr. Mate-pairs were --rf, but you don't see those much any more.
--rna-strandedness is a very different option, since it sets how reads are expected to align against genes. Most stranded protocols in use these days follow the dUTP-method, where read #2 in a pair has the same orientation as the transcript from which it arose. So either R or RF would typically be appropriate, unless the library is unstranded. In practice, I expect this is more useful if you plan to run stringTie downstream, since then the XS auxiliary tag is set appropriately.

I found this resource t visualise the relative orientation of reads. The --fr/--rf/--ff option is kind of disturbing... you say Mate-pairs were --rf, but you don't see those much any more. but In a "old" tool like tophat rf was not possible, when in a much recent tool they have introduced it...

There are some other sequencing types (i.e., not RNAseq) where different library types are now being produced. Presumably they wanted to allow some flexibility so hisat2 could be used with them as well.

source: Stack Exchange

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

sort: votes


text:
How to filter Ensembl cDNA and ncRNA FASTA files by primary assembly?
I'm currently performing differential expression analysis using alignment-free quantification using Kallisto. To do this, I need to create a Kallisto index using Ensembl's cDNA and ncRNA annotations available at the following two links:
ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/ncrna/Mus_musculus.GRCm38.ncrna.fa.gz
However, later on in my analysis I realized that I had multiple Ensembl gene IDs pointing to the same gene. For example, I was gettign results for both ENSMUSG00000026012 and ENSMUSG00000102412 which both map to Cd28. I realized that the reason for this error was due to the fact that the annotations I was using contained haplotypes.
I would like to remove these haplotypes, and to do this I believe that I need to rerun my analysis using a primary assembly. However, Ensembl doesn't have annotated cDNA and ncRNA FASTA files of the primary assembly. How can I filter these files to contain only genes present in the primary assembly? Is there a better way for me to solve this issue?

Both Ensembl and GenBank databases have their own ideas about what constitutes a gen, and the mapping of Ensembl genes to RefSeq genes is not one-to-one. There are some RefSeq genes that map to multiple Ensembl genes, some Ensembl genes that map to multiple RefSeq genes, and some genes that are not in the other database.

The simplest method is to download the GTF file for GRCm38 and filter that. You can then use one of the many tools out there (bedtools getfasta (make sure to enforce strandedness and possibly use the --split option, if you happen to use a BED12 file), gtf_to_fasta from tophat2 (it's the only reason to keep tophat2 installed), etc.) to make a transcriptomic fasta file.
Alternatively, Ensembl BioMart can be used to directly download a transcriptomic fasta file using a filter to select only those transcripts on pre-defined chromosomes. Other filters can be applied as well. If you don't already have a GTF file downloaded and one of the previously mentioned packages installed, then this may be the simplest way to get what you want.

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


text:
Program to make a haplotype network for a specific gene
Currently my lab manually uses pop art to make the network after using PHASE to id an individual's alleles. I would like to know how to generate the graph structure (i.e. the links between vertices), as well as how to draw the graph.
Input data are formatted as FASTA files with the sequence in nucleotide formatting.
Example:

I just feel like there's got to be a better, more time effective way to do this.

This is a haplotype map, where each node is proportional to the frequency of a given allele. The formal tool to do this is Network 5 found here. The problem with manaul haplotype maps is the area of the node should be proportional to the frequency, if this is not accurate it could be perceived to be misleading.
The approach for looking at each parsimonious state along a branch is a phylogenetics method which is performed via MacClade 4.0 ... requires an old Mac though.

What do network edges represent?

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


text:
Volcano plot in R
This question has also been asked on biostars
How can I reproduce this volcano plot?

I'm only able to do the traditional one, I'm kind knew too these field.

This plot is clearly done using core R functions. There are smoother alternatives how to make a pretty volcano plot (like ggplot with example here), but if you really wish to, here is my attempt to reproduce it :

I obviously had to generate data since I do not have the expression data from the figure, but the procedure will be about the same with the real data. Here is the R script :
### define interesting genes
interesting_genes <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
fc_interesting_genes <- c(-1.8, -1.3, -0.9, -0.85)
pval_interesting_genes <- log10(c(0.0005, 0.01, 1e-4, 4e-05))

### define colours (grey, blue, red)
palette <- c(rgb(0, 0, 0, max = 255, alpha = 65),
    rgb(0, 0, 255, max = 255, alpha = 125),
    rgb(255, 0, 0, max = 255, alpha = 125))

### generate some expression like data, fold change using normal distr, and tranformed uniform for pvalues
fold_changes <- c(rnorm(1000, 0, 2), fc_interesting_genes)
pvalues <- c(log10(runif(1000) / 3^(1 + abs(fold_changes[1:1000]))), pval_interesting_genes)

# find deferentially expressed genes to colour
left_biased <- fold_changes < -1
right_biased <- fold_changes > 1

point_type = 20
point_size = 1.5

### plot corpus
plot(fold_changes, pvalues,
     pch = point_type, cex = point_size, col = palette[1],
     xlab = "Fold Change (log2)", ylab = "P value",
     xlim = c(-6,6), ylim = c(-5, 0), yaxt="n")
aty <- axTicks(2)
labels <- sapply(aty,function(i)
            as.expression(bquote(10^ .(i)))
          )
labels[6] <- 1
axis(2,at=aty,labels=labels)

### plot deferentially expressed genes in colour + genes on interest with black
points(fold_changes[left_biased], pvalues[left_biased],
       pch = point_type, cex = point_size, col = palette[2])
points(fold_changes[right_biased], pvalues[right_biased],
       pch = point_type, cex = point_size, col = palette[3])
points(fc_interesting_genes, pval_interesting_genes, pch = point_type, cex = point_size)

### add text
labelx <- c(-4, -4, 4, 4)
labely <- c(-1.5, -0.5, -0.5, -1.5)
labels <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
text(labelx, labely, labels)
text(-4, -4, "down\n194 genes")
text(4, -4, "up\n156 genes")

### add lines between text and data points
for(i in 1:4){
    lines(c(labelx[i], fc_interesting_genes[i]), c(labely[i], pval_interesting_genes[i]))
}

Without a minimal reproducible example (MWE), I can't reproduce the plot but I would suggest using existing volcanoplot functions such as the limma package on Bioconductor. Here is a typical workflow for a differential expression analysis that produces a violinplot:
#install package
source("https://bioconductor.org/biocLite.R")
biocLite("limma")

#load package
library("limma")

#define phenotype/experimental conditions for samples
group = rep(c("Treatment","Control"),
c(number_of_treated_samples,number_of_controls))
design = model.matrix(~group);
colnames(design) = c("Control","TreatmentvsControl")

#fit empirical Bayes model
fit <- lmFit(NormalisedExpression, design)
fit <- eBayes(fit)

#summarise significant results
tt <- topTable(fit, coef="TreatmentvsControl",
                adjust="BH",n=nrow(NormalisedExpression))

#plot results
volcanoplot(fit, coef="TreatmentvsControl")

This plot can be customised in a similar manner to base R plots by passing the relevant arguments as shown in the limma documentation. This is the recommended plot format that readers in the field will be familiar with. While I do not recommend deviating from this standard, it is possible to reverse axes in R and the source code for the plot function can be retrieved from R by entering volcanoplot (without any arguments) if you wish to modify it.

Easy to visualize volcano plot using Python. Check this tutorial for inverted volcano plot

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


text:
Deseq2 output issue
3 Control replicate and 2 test sample ,so Im comparing my Control which is my Progenitor cells to Monocyte which is mature stage
this is my code im using
condition <- factor(c(rep("Control", 3),rep("Test", 2)),levels=c("Control", "Test"))

(coldata <- data.frame(row.names=colnames(countdata), condition))
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)

This is my deseq2 output
gene baseMean log2FoldChange lfcSE  stat pvalue padj CB1 CB2 CB3 GD5_1  GD5_2
ENSG00000144655 2030.6146124999 -5.0812043868   0.1183868766    -42.9203348727  0   0   3327.4890036209 3267.1877800456 3362.7407570783 98.5079578403   97.1475639145

Im not sure how come the pvalue and  padj is 0? even though the sample value there is a difference and there is also fold change ,something wrong with the result ?

What's your sample design? How many biological reps do you have?

I have updated it

A p-value of 0 is expected in this case. The p-values are calculated by comparing the Wald stat to a standard normal, so you can confirm that this is correct like this:
> pnorm(-42.9203348727)
[1] 0

It really means p < 2.225074e-308 See here for a more details

Yes i checked the post it has similar issue ,but pardon my ignorance but how come the deseq2 normalised value are so large CB1 CB2 CB3 GD5_1 GD5_2
3327.4890036209 3267.1877800456 3362.7407570783 98.5079578403 97.1475639145
yet these are not significant? can you explain me in simpler manner and relate with the biology , you can expand your current answer ..thank you in advance

small p-values are more significant, so a p-values less than 2.3e-308 is very, very significant. It's as significant as you can get


text:
What could be the reason for the samples not clustering?
I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in Stringtie manual. With the counts data I have used edgeR for Differential analysis.
Before that I have seen MDS plot showing in the following way:

I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.

Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?

Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.

Yes, you are right. If some thing went wrong what all should I check anything like mapped reads percentage or anything else?

Mapped percentages, size factors, and perhaps sparsity (percentage of genes with no counts). Basically you want to ensure that they're all similar. Also, try redoing the counts with featureCounts. Then you can eliminate an oddity in the counting as an issue.

Which is the mapping percentage in the above comments? It gives overall alignment rate.

The alignment rate is the mapping percentage.

Counts with featurecounts should I use bam or sorted.bam?

It doesn't matter which you use.

I initially checked on one bam file. Geneid, Chr, Start, End, Strand, Length, /path/to/A.sorted.bam. These are the columns in the featureCounts output. Is the last column read counts? Do you think this is fine? featureCounts -a annotation.gtf -F GTF -p -s 2 -T 8 -P -o countMatrix_all.txt /pathto/*.sorted.bam

Yes, the last column is the counts (I wish there were an option to simplify that output to just geneid and count). You probably don't need the -P option, not that it's likely hurting anything. BTW, that could be its own post here.

One last question. When I did samtools flagstat on one of the bam file i see this. 202671876 + 0 in total (QC-passed reads + QC-failed reads)
128036784 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
179886195 + 0 mapped (88.76% : N/A)
74635092 + 0 paired in sequencing
37317546 + 0 read1
37317546 + 0 read2
45092910 + 0 properly paired (60.42% : N/A)
48583902 + 0 with itself and mate mapped
3265509 + 0 singletons (4.38% : N/A)
621892 + 0 with mate mapped to a different chr
377778 + 0 with mate mapped to a different chr (mapQ>=5) In this what is mapped and properly paired?

Please post things like this as stand alone questions (make sure to google them first, this is answered elsewhere).

The largest variation differences in total expression will be seen in cell population differences. It's a good thing that your conditions aren't clustering on the first and second principal components when all genes are considered, because it means that cell population is probably not a factor in the variation of expression.
There may be some other principal component that does separate the conditions, but you might need to drill down fairly deep into the data. This will depend on how much the treatment conditions contribute to the total expression variation. We got lucky with our Th2 RNASeq study in that cell population was mostly on the first and second principal components, and treatment condition appeared on the third principal component.
If you want to force condition clustering, identify the genes that are differentially expressed, and only use those genes for the principal component analysis.

I'm not sure how useful doing a PCA on DE genes is. You already know your samples differ by group in those genes, after all.


text:
Regressing out unwanted sources of variation in single cell RNA-seq data
I have a dataset of single cell count data and I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes.
I know that count data is discrete data, and commonly follows a negative binomial distribution. To regress out these two confounders, should I use a linear model or GLM (poisson/negative binomial)? And how can I determine the optimal choice?

What do you mean by single cell count data? If you mean mRNA counts, then your data are compositional and should be treated as such.

yes mRNA counts, as I mentioned in the title. What regression technique is applicable for compositional data?

If I understand you correctly, "I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes" means that you want to mitigate the constrained geometry arising from multinomial sampling (aka the unit-sum problem). There are multiple transforms that can alleviate certain effects of compositionality (alr, clr, iqlr) or project the data into an unconstrained Euclidean space (ilr). The choice depends on subsequent analyses and their properties. ALDEx2, a widely used RNA-seq package, uses clr and iqlr. Take note, that they all assume certain properties.

Cleaning data before doing analysis can be more important than the method of data analysis. For UMIs, there's an obvious cleaning step that can be done: it would be better to filter out duplicates (prior to generating counts) than to try to incorporate that information in the count analysis.
I'm not familiar with single cell analysis, but I am aware that DESeq2 can be adapted to work with single cells, and incorporating additional variation is just a matter of adding an additional covariate into the DESeq2 formula model.

Thank you for your answer. In the meantime i found out that a linear regression is commonly used, but a negative binomial would be the favorable way to go.

Seurat has as part of its protocol a step where you filter based on UMI counts and percent mitochondrial
http://satijalab.org/seurat/pbmc3k_tutorial.html

Yes, but it uses a linear regression as a default. My question is specifically why a linear or a negative binomial regression.

If you just want to filter UMI related problems and cells with high mitochondrial to nuclear ratio, maybe there is no need to use modeling.

For UMIs: You can filter UMIs that are not supported by at least x reads (for instance 3 reads) and that probably arise from sequencing errors.
For high mitochondrial to nuclear ratio: You can compute (number of molecules of mitochondrial genes / number of molecules of nuclear genes) and plot this ratio for all the cells and try to put a threshold value for cells that have high values for this ratio. Thus, filter these cells that behave strange compared to the rest.

Thank you for your reply. I know these are standard threshold protocols, mostly dealing with badly sequenced cells or dying cells, however, this is not a correction for amplifcation bias. So I don't think this answer is what i'm looking for.

@DCZ I was following the Seurat tutorial for clustering cells (http://satijalab.org/seurat/pbmc3k_tutorial.html) when I found something very similar to your question. Literally, in the removing unwanted source of variation step: "In this simple example here for post-mitotic blood cells, we regress on the number of detected molecules per cell as well as the percentage mitochondrial gene content." It seems that they used a linear model, extracted from Buettner et al., 2015 (https://www.nature.com/articles/nbt.3102). Is this what you are looking for?. They even have an implemented function.

That's correct. The Seurat tutorial indeed does this, but the scaleData function uses a linear regression by default. As count data follows a negative binomial distribution, this seemed a peculiar choice to me.


text:
BioMart does not return a result
I am supplying to BioMart 81 ensemblIDs, but getting 80 results where 1 ensembleID: ENSMUSG00000083840 is absent. What does it mean if BioMart can not find the gene? Does it mean that it is a non-coding one or that there is just no symbol attached to it? Here is my BioMart query:
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","mgi_symbol","description"),values=genes,mart= mart)

Whenever you are wondering about things like this, just look up the identifier on the Ensembl web page. If you look up ENSMUSG00000083840, you will see:

This identifier is not in the current EnsEMBL database

In this case, it doesn't seem to have been replaced by anything and it wasn't a protein coding gene, so it looks like it was simply a wrong annotation which has been removed.


text:
Generate SNP/indel annotation in Arabidopsis
I have a genome-wide list of germline SNPs and short indels for Arabidopsis thaliana, which I generated with Varscan. Regardless of the tool used to generate them, I would like to annotate them, i.e. knowing which ones can cause an aminoacid change or an early stop, using the default Arabidopsis thaliana Columbia 0 cultivar, for which I have both the sequence (FASTA from TAIR10) and the updated annotation (GFF from Araport).
Chrom   Position    Ref Var
Chr1    626503  G   T
Chr1    926694  C   T
Chr1    5280350 C   A
Chr1    5699993 C   A
Chr1    7004559 G   A
Chr1    8325810 C   T
Chr1    9371723 T   G

What I want to do is similar to what Annovar does, but unfortunately Annovar does not support Arabidopsis. I was thinking of an already existing R pipeline that takes in a genome, an annotation, a SNP/indel list and boom, annotation. But I couldn't find any, except maybe snpEffect. Any tips?

Unfortunately I couldn't find a tool already out there.  I had to write something using python and sqlite.

The Ensembl Variant Effect Predictor will annotate your A. thaliana variants:
http://plants.ensembl.org/Oryza_sativa/Tools/VEP
The web based tool will use the reference genome assembly and gene annotation used in the Ensembl browser:
http://plants.ensembl.org/Arabidopsis_thaliana/Info/Annotation/#assembly
However, if you use the command-line version of VEP, you can provide your own reference genome assembly and gene annotation files to be used as a reference.
You can find all of the documentation here:
http://www.ensembl.org/info/index.html
Best wishes
Ben
Ensembl Helpdesk

I recently used SIFT4G to predict the deleteriousness for a SNP dataset in Arabidopsis.
I noticed that besides this, the tool also provides annotation for the SNPs (like gene region, amino acid substitution, start/stop etc.)
They have an Arabidopsis database readily available so you don't have to generate your own. Your input should be enough to run the tool, although you might have to do some reformatting (the tool requires VCF format files).
One thing though, not all SNPs are considered by the tool. Intronic SNPs and SNPs outside gene regions are ignored afaik.


text:
How are Principal Component analyses and Admixture analyses from a genetic alignment different?
How are Principal Component analyses and Admixture analyses from a genetic alignment different?
My understanding is that a PCA will take raw genetic differences across the entire alignment and plot them using dimensionality reduction techniques (single value decomposition for instance) whereas the Admixture analysis will plot likelihoods of sharing certain amounts of ancestry, based on allele frequency data. What I don't understand is how the two are fundamentally different. That is, isn't every variable position found in a PCA going to represent an allele used in Admixture?

The fundamental difference lies mostly in the math.
An admixture analysis involves assuming that genotypes (or more likely, genotype likelihoods) in an unknown sample can be modeled with Hardy-Weinburg equilibrium as arising from the combination from two or more source pools of background variants. The output is then an admixture proportion that dictates the mostly percent contribution from each of the source populations (e.g., 23% Irish, 15% Norwegian, 30% Northern African, ...). There are versions of this specifically for handling pecularities of NGS datasets.
A PCA is a very very generic method that uses the same input allelic information, but simply finds orthogonal projections of that onto lower dimensions that maximize the variance on the first principal component. You've probably seen plots like the following (from this paper about 1000 genomes projects data):

That allows you to get an idea about how populations are structured, but they're not so useful for determining how individual samples are most likely related to the various populations.

I find also a difference in the assumption of k ancestral population for admixture.
PCA/admixture differs in many aspects, but what about the difference between admixture and k-mean cluster - both likelihood based, assumption of k sources... the difference appears to be more subtle.  Hmm, maybe I should open a question.

Just to add to the great summary by Devon Ryan: admixture analysis tools are much more flexible than PCA (which is just a fixed mathematical operation), so they can be designed to incorporate LD patterns, phase information, different population evolution models etc.
Here's a paper from the creators of fineSTRUCTURE and Chromopainter admixture tools detailing the action of both - it seems they provide some comparisons with PCA as well:
http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002453


text:
Why chose CNN for a variant caller
Google released their variant caller DeepVariant which won the highest SNP performance award in the Precision FDA Truth challenge (99.999% accuracy).
From the linked github repo, we see that DeepVariant is a CNN, we provide images of the reads as an input and train the neural network with a stochastic gradient descent algorithm.
To me this is unintuitive. Why use images of reads when you already have the reads data?

https://en.wikipedia.org/wiki/Law_of_the_instrument

Where did you see they need images? The documentation says the input is a bam file, just like any other caller.

@terdon Images (of pileups, of course) are mentioned in the preprint. I think they're also using the Inception CNN, which is heavily used in image processing. I think many people forget that images are just matrices.

@terdon It seems that the first version was on images: We have made a number of improvements to the methodology as well. The biggest change was to move away from RGB-encoded (3-channel) pileup images and instead represent the aligned read data using a multi-channel tensor data layout. and now they use something else. Thank you!

As @DevonRyan has already said, these are not images per se. From a machine's point of view an image is just a bunch of matrices (one per colour channel). I'd say this is just an example of bad/confusing use of terminology.

It's more "images of reads" rather than actual images of reads. In reality, they're feeding in a pileup, which is just a matrix of letters or numbers (you can visualize this as an image, of course). This is more apparent if you look at one of their example python notebooks, CNNs are popular in the image processing field, but at the end of the day images are just large matrices of values so any time you can convert your input into that sort of format then a CNN can be useful (provided you have sufficient training data and some GPUs to throw at things).

I find the answer confusing: are you saying that the CNN is not fed images? Because the way google describes DeepVariant, they are indeed creating pileup image data and putting that through their CNN.

It's not an image, it's an "image". The DeepVariant are a bit ambiguous on that point, though at the end of the day it's a multidimensional matrix of values (the relation of that to an image being obvious).

Why use images of reads when you already have the reads data?

Neural networks take one or multiple N-dimension arrays (n-d arrays) as input; they can't take arbitrary forms of data. To feed reads and their alignment into the network, you have to encode them as n-d arrays. If you ignore "colors" for now, a straightforward choice is a (windowLength, readDepth) 2D-array as you see in an alignment browser – it is an image. The "colors" of this image are often called "channels" or "features" in deep learning (DL). They represent base types and base quality in deepVariant. In the end, you get a (windowLength, readDepth, nFeatures) 3D array as input. This is the most natural encoding to DL researchers.

Why chose CNN for a variant caller

That is because CNN has been repeatedly proved to be the best architecture for such (height, width, nChannels) 3D arrays as input. Read alignments are exactly that.

On a historical note, deepVariant was first built on top of google's now-deprecated DisBelief framework that only takes RGB images as input. They had to squeeze all types of information into three color channels. Now with TensorFlow, an image can have an arbitrary number of color channels. This simplifies implementation and improves accuracy. In addition to this image, deepVariant also takes other information as input.


text:
Mugsy error: Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501
I have installed mugsy in order to create a multiple genome alignment and a phylogenetic tree of several species of nematodes. The following command successfully pulls out the help:
output/software/mugsy_x86-64-v1r2.3/mugsy -h

To test that the software is installed correctly and that I am using it correctly, I ran the following alignment between the Wormbase sequence and the Ensembl sequence of the C. elegans genome:
output/software/mugsy_x86-64-v1r2.3/mugsy  --directory /.../project/output/software/mugsy_x86-64-v1r2.3/ --prefix mygenomes output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS263.genomic.fa output/genome/celegans/ref/ensembl/seq/celegans.fa

I get the following output:
Parsing sequences for c_elegans  num_seqs:7
Parsing sequences for celegans  num_seqs:1
Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501.

The headers of the wormbase sequences are:
>I
>II
>III
>IV
>V
>X
>MtDNA

and the headers of the ensembl sequences are:
>I dna:chromosome chromosome:WBcel235:I:1:15072434:1 REF
>II dna:chromosome chromosome:WBcel235:II:1:15279421:1 REF
>III dna:chromosome chromosome:WBcel235:III:1:13783801:1 REF
>IV dna:chromosome chromosome:WBcel235:IV:1:17493829:1 REF
>MtDNA dna:chromosome chromosome:WBcel235:MtDNA:1:13794:1 REF
>V dna:chromosome chromosome:WBcel235:V:1:20924180:1 REF
>X dna:chromosome chromosome:WBcel235:X:1:17718942:1 REF

Maybe this line:
    Parsing sequences for celegans  num_seqs:1 should make me worry. Why does mugsy think there is only one sequence in that genome? Why do I get the following error and how can I fix it?:
Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501.

I am working on a GNU/Linux platform.


text:
Why models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decrease
I am reading Elowitz et al. (2002). When talking about intrinsic noise, i.e. noise due to microscopic events that govern which reactions occur and in what order, it is stated that:

Models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decreases.

Noise, defined in the article as the standard deviation divided by the mean, can be divided in intrinsic (due to microscopic events during transcription) and extrinsic (due to localization, state and concentration of other cellular components, e.g. polymerases).
I think the intrinsic noise should increase at low level of transcript amount because the average (mean) amount will be lower but deviation decrease is smaller than mean decrease. I don't know, however, if this is the actual reason (and why deviation decrease should be smaller than mean decrease) or if deviation decreases with the mean and there's another more likely explanation.

Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S. Stochastic gene expression in a single cell. Science 297, 1183–1186 (2002).

Reading Kaufmann and van Oudenaarden (2007), it seems to validate the first alternative (using results from the Central Limit Theorem):

Although biochemical fluctuations influence all stages of gene expression, those involving molecules in extremely low abundance are expected from a statistical standpoint to be larger in magnitude and therefore to contribute disproportionately to the overall variation between cells.

Kaufmann, B. B. & van Oudenaarden, A. Stochastic gene expression: from single molecules to the proteome. Curr. Opin. Genet. Dev. 17, 107–112 (2007).

The inverse relationship between variance and mean is a consequence of the stochastic nature of gene expression, which can be modeled as a discrete birth-death process. One of the simplest models assumes that the gene is expressed constitutively. Then, using the chemical master equation, mRNA production can be shown to evolve as a discrete Poisson process. The steady-state distribution is Poisson. For the Poisson distribution, coefficient of variation ($\sigma^2/\mu^2$) is equal to $1/\mu$. Further, if protein is translated in bursts from mRNA (translation burst), protein noise can be written as
\begin{equation}\label{eq1}\tag{1}
\eta^2_P = \frac{\sigma_P^2}{\mu_P^2} = \frac{B}{\mu_P},
\end{equation}
where $\eta_P^2, \sigma_P^2,\mu_P,$ and B are the coefficient of variation, variance, mean and burst size for protein. From \eqref{eq1}, the inverse relationship is evident.
Gene expression models can be made more detailed by assuming different states for the promoter. The promoter switches between these states, and each of these states can have different rates of transcription. For these models, $\eta_P^2$ can again be shown to follow an inverse relationship with $\mu_P^2$.
Reference:
Paulsson, Johan. "Models of stochastic gene expression." Physics of life reviews 2.2 (2005): 157-175.


text:
Obtaining identifier from plot of sequence GC%
I was looking at this script from https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297 and wanted to use it to determine outliers for tentative Horizontal gene transfer.
The data set is provided in the example, under the section 20.2.2 Plot of sequence GC%.
I want to yield a identifier for the x axis instead of just an arbitrary number of genes or cds' so I can actually determine which sequences have the outlying values, as I can not with the current graph.

I am not great at python or understanding list comprehensions, so if anyone could augment the code provided in the example to do so, it would be appreciated.
I provided the code below also.
#from : https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297, data set under
#section 20.2.2 Plot of sequence GC%, ls_orchid.fasta

from Bio import SeqIO
from Bio.SeqUtils import GC

gc_values = sorted(GC(rec.seq) for rec in SeqIO.parse("ls_orchid.fasta", "fasta"))
import pylab
pylab.plot(gc_values)
pylab.title("%i orchid sequences\nGC%% %0.1f to %0.1f" \
            % (len(gc_values),min(gc_values),max(gc_values)))
pylab.xlabel("Genes")
pylab.ylabel("GC%")
pylab.show()

If you goal is to guess putative HGT, is there a reason why not to use one of the pipelines built for it?

I was not aware of the pipelines. Or i guess  i didnt feel comfortable using them. @KamilSJaron      did you have any suggestions?

The tricky bit here is that you need a list of sequence records sorted in the same order as the CG values. You can accomplish this by making tuples of the seq records and GC values, then sorting the tuples specifying the GC values as the sorting key:
gc_w_records = sorted(((GC(rec.seq), rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")), key=lambda tup: tup[0])
gc_values = [x[0] for x in gc_w_records]
records = [x[1] for x in gc_w_records]

I don't know much about matplotlib and it's going to be challenging to display 100 IDs, but you'll probably be interested in the extremes of the GC distribution, which you can obtain directly from the sorted list of records:
records[0:2]
# [SeqRecord(seq=Seq('CGTCACGAGGTCTCCGGATGTGACCCTGCGGAAGGATCATTGTTGAGATCACAT...CAT', SingleLetterAlphabet()), id='gi|2765587|emb|Z78462.1|PSZ78462', name='gi|2765587|emb|Z78462.1|PSZ78462', description='gi|2765587|emb|Z78462.1|PSZ78462 P.sukhakulii 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...AGG', SingleLetterAlphabet()), id='gi|2765568|emb|Z78443.1|PLZ78443', name='gi|2765568|emb|Z78443.1|PLZ78443', description='gi|2765568|emb|Z78443.1|PLZ78443 P.lawrenceanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

or:
records[-3:-1]
# [SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet()), id='gi|2765656|emb|Z78531.1|CFZ78531', name='gi|2765656|emb|Z78531.1|CFZ78531', description='gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGCTCATTGTTGAGACCGCAA...AAG', SingleLetterAlphabet()), id='gi|2765625|emb|Z78500.1|PWZ78500', name='gi|2765625|emb|Z78500.1|PWZ78500', description='gi|2765625|emb|Z78500.1|PWZ78500 P.warszewiczianum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

+1 but would be marginally faster to use operator.itemgetter(0) as the key for sorted, or just omit the key and the sorting will occur by first the GC then the record

I did try that, but you can't sort sequence records: "NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest."

Ah yes good point, forgot about that, itemgetter(0) would be the way to go then, although by not essential

There are various attempts to guess putative horizontal gene transfer using sequence composition (kmer-frequencies or GC content), you do not have to implement a new method by your own.
I co-authored one of them called SigHunt. It's kind of working, but the code is really messy and there are now bit more precise methods out there. The advantage is SigHunt is that it provides access to all the data calculated on the way (the tetranucleotide frequencies), so if you want to play with it, it's fun.
The method I would probably go for now if I would have searched for HGT is MTGIpick. The most of the tools detects HGT within a genome by a sliding window, but MTGIpick is refining boundaries where the change of kmer frequencies is rapid. I never tried it personally but the approach is sound to me.
Alternatively there are other paradigms for detection of HGT. Phylogeny based methods (there are many, but for instance) or by mapping reads to species boundaries. Phylogeny based methods usually provide stronger evident, but have a disadvantage to be database dependent (i.e. you can not build a tree if you don't have sufficient coverage of homologs across tree of life).
I would also like to mention that right now I do not follow the HGT literature closely, there might be some big new I just missed.


text:
Is there a publicly available tumor-normal sample?
I am looking for a publicly available matched tumor-normal sample. I need Illumina fastq reads (or an aligned bam file, since I could extract the reads from it) from a tumor and a matching, non-tumor control set from the same patient. Ideally, I would also have a truth set of known variants in the sample, but I can live without that.
I checked the Next Generation Cancer Knowledge Network but although they do seem to have appropriate bam files, which I could convert to fastq, all of the ones I found were "controlled access" and getting access seems very complicated.
Is there something like the public NA12878 sample for tumor-normal analyses?

The TCGA does provide the raw fastq too, and for some patients one has sample from the tumor and from a healthy region. But I am not sure it is public as without registration

@Llopis as far as I can tell, the registration process is incredibly convoluted and involves getting the "Signing officer" of the organization (PI, CEO, whatever) to register the organization with two separate governmental bodies (eRA and then dbGAP) and then ask for access. It might be possible for me to ask my boss to do this, but I would much prefer to find more easily accessible data.

Most patient data is stored in dbGAP because of privacy issues, so you'll struggle to find anything that's unrestricted.

@heathobrien yeah, which is why I've been struggling :) I was hoping there was a Henrietta Lacks analogue for the modern day or that NIST had released a tumor-normal pair or something.

I thought the HeLa genome got pulled down for exactly that reason. Is it available again?

Not that I know of, no. Sorry, I was thinking of her cells and how they've been used for so many years, not about her genome. I come from the era before NGS.

If you work at a university, getting access to dbGap protected data is no where near as difficult as it seems. Getting registered on the system originally took about an hour of my time and about half a week of waiting. When I apply for a new dataset it takes me about half an hour per project.

@IanSudbery unfortunately, I am no longer in academia, so I would need to ask the CEO of the small startup I work at to do this, and he barely has enough time to breathe as it is.

Please take a look at
An open access pilot freely sharing cancer genomic data from participants in Texas
Although I have not worked on this paper but author claimed that both tumor and normal data is publicly available. In case if you still need something more, than you can work on virtual normal data. And there are several papers which provide machnine learning techniques to overcome the deficiency of normal data (not much efficient but still they give considerable results). Hope this information will help you. Thanks

Argh, this looks perfect, but the site seems to be down! I have emailed the author of the paper to see if it's just a temporary glitch. Thanks!

The site is no longer maintained, but I contacted the authors and they very helpfully put it back up for me!

@terdon: Then you should look for virtual normal database than or, if interested, machine learning techniques.


text:
Count files using htseq-count?
I am doing RNA-Seq data analysis using data available on SRA (SRR6047326). Samples are from the pig and they are paired. So, I did alignment using HISAT2.I used pig reference genome available on ENSEMBL. Then, I sorted .sam output files obtained from HISAT2. Now I am trying to get the count file using sorted sam file and it is not working. I am using annotated GTF file for pig which I downloaded from Ensembl.
I get count files using following commands, but all the counts are 0.
htseq-counts -i no myfile.sam genes.gtf>output.count

What is wrong here ? Are there any alternatives to htseq-count?

Help us to help you. What does "not working" mean? Please post relevant error messages, if any.

I runs well and gives me the count file. But there are zero counts for all the genes. At the end of count files it scores most of the genes are on --alignment_not_unique.

I tried kallisto and it seems to work fine with kallisto.I can see counts for the multiple genes in tsv file generated by kallisto.

Is that the actual command you ran? Because the "no" in -i no myfile.sam is odd, and since kallisto is returning some reasonable counts we can assume the issue with the htseq command and not the alignments (which was my 1st thought).

Speaking of kallisto, it is also fine to use those results for your analysis, unless that there a specific reason to use htseq-count.

And I gave a shot to "Salmon" too. It is awesome. Basically it works like Kallisto. Thanks for your help guys.

I m not sure about htseq-counts but I would suggest go for featureCounts, I will update my answer and share the script/command. Give it a try
the command i use
featureCounts -T 50 -C -p -t exon -g gene_id -a /run/media/data2/gencode.v21.annotation.gtf  -o PRIMARY_CELL.txt $(ls *.bam | sort )

Thanks. I haven't use featureCounts. I found that it could run on R. Hopefully it works.

yes you can but it takes lot of memory

I am running on cluster. But, I have a question here. I have .sam files obtained from HISAT2 which were generated from the paired fastq files. Can I use that .sam files directly for feature counts ? Do I need to sort the sam file and then use feature counts ?

I used the following commands
fcounts <- featureCounts(files=“reads.sam”,genome="genes.gtf").
And my results were:
Unassigned_NoFeatures = 5305153, and Unassigned_Unmapped=49253221. So, none of the genes were mapped. Any idea what could have gone wrong here ?

i have updated my answer see if that helps becuase i never use the sub-read package in R as such , i think you need to convert it into bam files...that would be better...

For some reason I was not able to get any of my reads mapped. But I was able to do my analysis using Salmon.

it shouldn't be case may be some minor issue because i use star to align data and then to get reads i use featureCounts ,in star i get bam sorted file as my output but I don;t think if you don;t get sorted files you wont be able to get reads using featureCounts ..but anyways it worked may be you can give it a try again..

Yes, I found a lot of papers using STAR aligner. I will give it a try. Thanks.


text:
well identifiers in a 1536-well mtp plate
There's a commonly accepted convention for naming the wells in a 384-well-plate:
A01..A24
B01..B24
..
P01..P24

It's concise and unambiguous and thus very useful.
For 1536-well plates (32x48) it's more difficult because there's more rows than letters in the alphabet. It seems there's no standard.
Which nomenclatures are out in the field and which is the best?

could someone with more reputation than me create and add the tags screening and mtp? and maybe modify the tag hts so that it includes high-throughput-screening?

Can you perhaps modify this such that it's on-topic here? As is, it seems neither directly related to bioinformatics nor a question with a best answer (the goal of the stackexchange network).

the goad is to find the most suitable scheme for naming well ids in a piece of bioinformatics software, so I don't really understand your concern. Is high-throughput-screening not a part of bioinformatics?

Screening is relevant, but random "how do I label wells" questions aren't necessarily.

Looking into standards/conventions of how to organize biological data most certainly is relevant to bioinformatics, and best answers can be envisioned: "There is now standard X", "Commonly used convention Y", "There is none yet, but Z has proven useful". I'm glad this helpful question was not closed/deleted.

Here's one way I found:
A01..A48
..
Z01..Z48
AA01..AA48
..
AF01..AF48

Pros:

Excel does it like that with columns
it's clear what's the row and what's the column

Cons:

alphabetic sorting of well ids does not what you want
when reading A01 it's not clear whether the plate size is 384 or 1536

and here's a scheme I made up myself:
AA01..AA48
..
AZ01..AZ48
BA01..BA48
..
BF01..BF48

Pros:

alphabetic sorting of well ids does what you expect
it's clear what's the row and what's the column
it's easy to distinguish well ids of 384-well plate from those of a 1536-well plate

Cons:

it's not used anywhere AFAIK


text:
Script to run everything in a loop for extracting tar.gz files into fastq and to bam with alignment?
Original tar.gz files look like below:
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz

Before extracting the tar.gz files, I changed their names based on names.text file as below and using the command:
xargs -n2 mv -i <<<"$(sed '1d; s/$/.tar.gz/' names.txt)"

Tar.gz_filenames       Sample_Names
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz                Sample2
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz                Sample1
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz                Sample3
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz                       Sample4

After changing the names, looks like below:
Sample1.tar.gz
Sample2.tar.gz
Sample3.tar.gz
Sample4.tar.gz

Using the tar.gz files after changing names I used the below code:
for i in /doc/*.tar.gz
do
    tar xvzf $i -C $TMPDIR
    for sample in $TMPDIR/*1.fastq
    do
        dir2="/destination/directory"
        base=$(basename $sample ".1.fastq")
        module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq | samtools view -Sb - > ${dir2}/${base}.bam
    done
done

When I run the above code, the output is like below:
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_2.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_2.fastq
Warning: Could not open read file "/scratch/username/slurm-job.36320423/140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1

I dont mins whatever the name of fastq's in $TMPDIR but from $TMPDIR hisat2 should use those fastq's and give output bam files in /destination like below:
Sample1.bam
Sample2.bam
Sample3.bam
Sample4.bam

I would like to run a script doing tar.gz files -> fastq -> bam (everything in loop)

This should work for you:
for i in *.tar.gz
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

(If this is the sort of thing you'll be doing often, it would be worth investing the time to learn a more robust approach such as Snakemake or Nextflow).
Edit:
If you want to run on all files at the same time, the easiest way is find / xargs
Change you script to this:
for i in $@
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

Then run it like this:
find . -d 1 -name '*.tar.gz' | xargs -n 1 bash NAMEOFSCRIPT

This works for SGE. I'm not sure about the details of SLURM. You might have to figure out how to run each job in the background or use gnu parallel or something. As @Ian Sudbury says, pipeline software can take care of all these details for you.

Definately +1 for putting all this in some sort of pipeline. Many of these  will now handle automated distribution of many files to your SLURM cluster. Job tracking. Restarts after errors etc.

@heathobrien Hi, The above code uses one after another tar.gz file. It is taking huge time. How should I make the code work for all the tar.gz files at a time?

see updated answer

Even I haven't worked on SLURM before. Do you think if I give #SBATCH --array=1-6%3 will work? lets say I have 6 tar.gz files so 3 jobs at a time.

Your samples end with _1.fastq and _2.fastq, so you need to change this:
-1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq

to this:
-1 $TMPDIR/${base}_1.fastq -2 $TMPDIR/${base}_2.fastq

Thanq. ya this was the only error in my code.

Hi, The above code uses one after another tar.gz file. It is taking huge time. How should I make the code work for all the tar.gz files at a time?

Untar the files and use a workflow manager like snakemake. BTW, make sure you have sufficient resources to process many samples in parallel.

But as soon as Untar done for one file it should be the input for hisat2. Yes, I'm using 20 cpu cores, 1 week time and tmp 500G, memory usage 4G.

Then put the untar command into something like Snakemake too.

Never used a snakemake till now. May I know how I can give that?

Please google such things and post a new question if you have further questions.


text:
running a samtools command for multiple bam files from 1000 genomes project
I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to download all of the bam files for the entire genomic region, just those intervals which I am interested in them.
I guess I can use this (I found from biostars forum)
samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)

However, I am just a bit confused. First of all, when I try the code above directly from the command line, I get an error message "No such file or directory"! Am I doing something wrong? or should I provide a specific path?
Second, I want to run  it for the mapped files for all individuals: *.bam”. through a loop.
Can anyone help me to solve my problem?

The command you're looking for is:
samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

For a series of intervals, you'd be better off scripting something with pysam:
import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

You might see if there are bigWig files available, they'd be faster to process.

,thanks. I could rum the sam tools command for a single file. Thanks for the pysam, I will try it. But do you know how can I run a samtools command through a loop? just to remind that I do not have bam files downloaded.

You can do for loops in bash. It can be pretty painful to script anything complicated in bash, though.

yes, I know the for loop solution. my problem is that the bam files for each individual are in separate folders. I am trying this: for file in http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data//alignment/.bam;
do /data/programs/samtools-1.3.1/samtools view -c "${file}" 2:1000000-2000000
done       but get an error message.

You can't wild-card URLs, you'll need to provide them as a list.

@Anna1364 please add all this information to your question. We can't answer you if you don't tell us all the details. Give us a list of the files you want to read and explain what you need to do with those files.


text:
Select top 100 genes ranked by variance in read counts
I have a gene expression data (read counts) table with many samples. Showing some of them below:
Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220

I want to select top 100 genes with reverse sort of variance. How to do that?

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it m:
sel = order(apply(m, 1, var), decreasing=TRUE)[1:100]

sel then contains the indices into your matrix (or the original dataframe).
BTW, I hope you're not filtering by variance before performing differential expression. That would be incredibly incorrect. Doing this sort of filtering for performing PCA or something like that can be useful, though you should really be using normalized counts (or a transformation of them).

Thanks. yes, I'm doing this filtering for PCA not for differential expression.

BTW, before doing this I converted Ensembl_ID to gene names and I have multiple genes now so I did like this U2 <- aggregate(. ~ gene_name, data = U, max) and on that I will be calculating variance. Do you think this is right?

Usually what I do is using the mean of the different transcripts for each gene, but it is right (if you can justify it).

@raju Do not convert to gene names until you are completely finished with the analysis.

Ok. But if there are duplicate Ensembl_ids can I do this U2 <- aggregate(. ~ Ensembl_ID, data = U, max)

And when I did this sel = order(apply(m, 1, var), decreasing=TRUE)[1:100] I see int[1:100] with numbers.

If there are duplicate ensembl IDs then you did something wrong. You should see a vector with 100 numbers, as I said their the indices into m.


text:
Help Adjusting Bedtools Output in R
I have been asked to write a script in R to edit, using data.frame, this .fa file that is the result from extracting gene sequences from a repeat masker coordinate file on the genome I am working with. Each input looks as such:
>rnd-5_family-5445_Unspecified(-)
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_Unspecified(-)
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATA
AAAATCTTTCTACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTG
TAATTTACAAAATGAATTACTTGAAGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCT
ACCGATGCCAAAAATTGACGATATTCTATTATTTAATCTAAATGTTAGTCTTTACATACCCTCCCCTA
ATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGCTAAATCTAATTCAATAGATT
GATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATTGCATATAA
CTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATG
AAAGTGATC
>rnd-5_family-2614_Unspecified(-)
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATT
AGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCT
ACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAA
ATGAATTACTT
>rnd-4_family-798_Unspecified(-)
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAA
ATAAACCTTATACTGCTAT

This is the repeat masker coordinates that the above were derived from:
>KB824701.1 417 478 rnd-5_family-5445_Unspecified . -
>KB824701.1 587 1072 rnd-5_family-2614_Unspecified . -
>KB824701.1 914 1129 rnd-5_family-2614_Unspecified . -
>KB824701.1 1138 1225 rnd-4_family-798_Unspecified . -

and ideally I would like it to be in the format of this, for example using the first one:
>rnd-5_family-5445_KB824701.1_417_478

I was told to use data.frame in R but I am not sure exactly how this will work? Can someone help please

Hi Noob, what have you done? Are you familiar with R or other shell programs? This seems as homework could you confirm it is so? I wouldn't like to post a solution to get you later accused to plagiarism... or worse, to provide an answer and you just used it without learning from it.

It is not homework, it is actually a project I am being to asked to help work on while also learning to skills so I have a lot of trouble trying to work and learn at the same time. I used a bedtools command getfasta to extract the sequences from a genome I was given repeat masker coordinates for already

Just checking that I understood correctly, you want to change from >KB824701.1 417 478 rnd-5_family-5445_Unspecified . - to >rnd-5_family-5445_KB824701.1_417_478? Or from  >rnd-5_family-5445_Unspecified(-) to >rnd-5_family-5445_KB824701.1_417_478. Are you comfortable with a other solutions not involving R ?

In case an answer would still be relevant after a year and a half, below is a solution with the Biostrings package.
library(Biostrings)

# loading data.table for the fread() function.
# fread() is not only fast but is quite good at guessing how
# the file is formatted
library(data.table)

# readDNAStringSet() function creates a "DNAStringSet" object
my_fasta <- readDNAStringSet("fasta.fa")

> my_fasta
  A DNAStringSet instance of length 4
    width seq                                                                names
[1]    61 TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG      rnd-5_family-5445...
[2]   485 AATCTGAGATTAGATTTATTTATCTGTTATCC...GGCCAATCATTCAAAAGTTATGAAAGTGATC rnd-5_family-2614...
[3]   215 TTTCTAGTTTGTATTTGAGGAAAGAAAATCCT...TAGAAATTGTAATTTACAAAATGAATTACTT rnd-5_family-2614...
[4]    87 ATGCTTAATTTGCTAAGCATTTAATAACCTTT...GTAGACATTAAAATAAACCTTATACTGCTAT rnd-4_family-798_...

# the fasta header is stored at the @NAMES slot of the @ranges slot of the
# DNAStringSet object created above
> my_fasta@ranges@NAMES
[1] "rnd-5_family-5445_Unspecified(-)" "rnd-5_family-2614_Unspecified(-)"
[3] "rnd-5_family-2614_Unspecified(-)" "rnd-4_family-798_Unspecified(-)"

# removing the trailing "_Unpecified..."
my_fasta@ranges@NAMES <- sub("_Unspecified.*",
                              "",
                              my_fasta@ranges@NAMES)

my_coord <- fread("repeat_masker.txt")
> my_coord
           V1   V2   V3                            V4 V5 V6
1: KB824701.1  417  478 rnd-5_family-5445_Unspecified  .  -
2: KB824701.1  587 1072 rnd-5_family-2614_Unspecified  .  -
3: KB824701.1  914 1129 rnd-5_family-2614_Unspecified  .  -
4: KB824701.1 1138 1225  rnd-4_family-798_Unspecified  .  -
#removing the ">" sign

my_coord$V1 <- sub(">",
                   "",
                   my_coord$V1)

# appending first 3 columns of the coordinate file to "NAMES" as requested by the OP
my_fasta@ranges@NAMES <- paste(my_fasta@ranges@NAMES,
                               my_coord$V1,
                               my_coord$V2,
                               my_coord$V3,
                               sep = "_")

# writing the modified DNAStringSet object into a fasta file
writeXStringSet(my_fasta, "edited_fasta.fa")

>rnd-5_family-5445_KB824701.1_417_478
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_KB824701.1_587_1072
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTA
CCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTTGA
AGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCTACCGATGCCAAAAATTGACGATATTCTATTATTTAA
TCTAAATGTTAGTCTTTACATACCCTCCCCTAATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGC
TAAATCTAATTCAATAGATTGATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATT
GCATATAACTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATGAAAG
TGATC
>rnd-5_family-2614_KB824701.1_914_1129
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATTAGATTTATTTAT
CTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTACCACTTTTGATTCATTTTTTAAA
TCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTT
>rnd-4_family-798_KB824701.1_1138_1225
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAAATAAACCTTATA
CTGCTAT


text:
SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why?
The SAM file format specification makes clear that optional fields must follow the format TAG:TYPE:VALUE
Type   Regexp matching                                       VALUE Description
A     [!-~]                                                  Printable character
i     [-+]?[0-9]+                                            Signed integer
f     [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?                 Single-precision floating number
Z     [ !-~]*                                                Printable string, including space
H     ([0-9A-F][0-9A-F])*                                    Byte array in the Hex format10
B     [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+    Integer or numeric array

However, looking at the samtools source code, I noticed this:
https://github.com/samtools/samtools/blob/develop/sam_view.c#L920-L931
switch (aux_type) {
    case 'H':
    case 'Z':
        if (kputs(bam_aux2Z(s), linebuf) < 0) return false;
        break;
    case 'i': kputw(bam_aux2i(s), linebuf); break;
    case 'I': kputuw(bam_aux2i(s), linebuf); break;
    case 'A': kputc(bam_aux2A(s), linebuf); break;
    case 'f': kputd(bam_aux2f(s), linebuf); break;
    case 'B': kputs("*** Unhandled aux type ***", linebuf); return false;
    default:  kputs("*** Unknown aux type ***", linebuf); return false;
}

Questions:
(1) I have never seen an "Integer or numeric array" of case "B". What does such an tag look like? Is there an example?
(2) Why is this currently "*** Unhandled aux type ***"? Was this the case with previous versions of samtools?

An example would be ZC:B:i,271,271,1,0, which is from TMAP. The i part denotes the values held in the array, which are comma-separated.
That code is only used in samtools fastq, where that tag type isn't output into fastq files.

Thanks for the help. "which is from TMAP"---where did you find this?

Also, how command is this tag? I don't think I've seen this type of tag "in the wild" in any BAM files....though perhaps I'm very wrong.

That particular example can be found here. You don't see tags like this that frequently, I had to search for a concrete example since I don't recall seeing any locally in a few years.

Last questions: (1) I haven't seen this either. Is there any evidence that this tag format has largely been abandoned in recent years? (2) Is there a reason why the tag isn't output into FASTQ files?

That's used under the hood if you have more than 65535 CIGAR operations, but otherwise hasn't really seen wide-spread use. I guess no one ever saw the point in outputting that in fastq files. Presumably someone could make a PR to add it.

"if you have more than 65535 CIGAR operations" Could you detail this a bit more?

I'm not sure what there is to detail there. Really long noisy reads can result in a lot of CIGAR operations, which exceed the bounds of the actual CIGAR field in a BAM file. Thus the CG auxiliary tag, which is also an array at its core.

Thank you. I just needed a further clarification. I think I understand now, thank you


text:
How can I find the chromosomal location of a list of genes?
I have a list of genes nearly 20000:
gene name (column1) and coordinates (columns 2 and 3)
head genes
    Ma0_73Ns    2359    2923
    Ma0_73Ns    4727    4868

I want to find the chromosomal location of each gene, on which chromosome each gene located? my study organism is Helianthus annuus (sunflower)

From what organism?

Do you know what the coordinates you mention refer to?

Biomart is the best tool for that. You don't even need any script for that. The user interphase is pretty neat and easy. You can refer to this tutorial for easy understanding - Click to watch
Regards
Arman

Helianthus annuus is an Ensembl plants genome. You should be able to use the bioMart there to pull out genomic co-ordinates (chromosomes etc) for genes. However, searching for the gene names you posted does not return any results. The same is true of the two websites posted above. Are they real gene names or just examples you made up for the post?

Ok so I found two websites that could help. The one I found first is more intuitive to go through, but it looks like it could be outdated as the second website I found seems more recent so it probably has a more complete genome to use as the reference. (Sorry, I haven't worked personally with sunflowers and don't know for sure where you would go for the "official" reference genome, or if there even is an "official" one, but hopefully one of these websites will help.)
1) This is a website I found that looks like it is a website for the sunflower genome and includes gene annotation information. (The recent publications list only goes up till 2015 though so it may not be the most recent data)
https://www.sunflowergenome.org/annotations-data/
If you scroll down to the bottom of this page there is a text file detailing the descriptions of how the files are formatted.
For example, the linked to file gives the following:

Ha412v1r1_CDS_v1.0.fasta.gz

Description: Nucleotide CDS for each gene
Header format: >TranscriptID GeneID Chromosome_start_end_strand
Example: >Ha1_00044280-RA Ha1_00044280 Ha1_94590914-94591860_-

If you download this fasta file you could then create a python script that will read in the file and extract the headers. If each header is a string variable called "mytext", you could then extract the TranscriptID, GeneID, and chromosomal info and using mytext.split(" "). You could then split on "_" to further parse the chromosomal location info.  You could save the extracted information in a dictionary using the gene name as a key and parse through the original file you mentioned containing the list of genes, in order to get the chromosomal location information you desire.  The one caveat is that this is assuming that your gene names are the same as the geneIDs the above website uses.
I did not exhaustively examine the list of files and perhaps there is another one that would be a better one to parse through for your purpose, but the linked to text file gives more details on the files and their formats and presumably they could be parsed in a similar manner.
2) Double checking my first answer I found this article:

Badouin H, Gouzy J, Grassa CJ, Murat F, Staton SE, Cottret L,
  Lelandais-Brière C, Owens GL, Carrère S, Mayjonade B, Legrand L. The
  sunflower genome provides insights into oil metabolism, flowering and
  Asterid evolution. Nature. 2017 Jun;546(7656):148.

Which links to a website as follows:

"The integrative web interface Heliagene provides visualization, querying tools for data mining and network exploration for the community (https://www.heliagene.org)."

So it looks like this website also has reference genomes for sunflowers, although a lot of the annotation information seems to be locked unless you have a login. Like I mentioned previously, this website is a little less intuitive to go through but you might be able to find something helpful here if you think the other website isn't official or recent enough.

If you can add a column of RefSeq identifiers for each gene to your table, first download a whole-genome GFF file for your organism:
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/4232/101/GCF_002127325.2_HanXRQr2.0-SUNRISE/GCF_002127325.2_HanXRQr2.0-SUNRISE_genomic.gff.gz

and then open the file in Python using a package like gffutils and extract all the rows which contain a RefSeq identifier located in table (the RefSeq identifier is the first column in the gff file linked in the wget command above). For all the rows you extract, which chromosome it belongs to will be stated in the Attributes column, which you can also extract using gffutils.


text:
Official documentation for textual BLAST output?
This is pretty vexing: I got asked to explain the meaning of the match line in a blastp output, especially the subtle difference between a blank and plus. E.g.:
 Query  181  NHIEEL...
             N I+EL...
 Sbjct  181  NQIKEL...

For me this is under the "everyone knows this" category, but I wanted to have an authoritative reference I could cite.
Alas, none of the documents on the BLAST help page seems to have this information. My Google fu also failed as I was not able to find an authoritative source which I could cite.
Also, I'm not sure whether BLAST introduced this notation or it simply picked it up somewhere else.
Any pointers appreciated.

The reference that is mentioned when I run a BLASTP search using NCBI's web service is this one. Here's some relevant text from the caption of Figure 3 of that paper:

On the central lines of the alignment, identities are echoed and substitutions to which the BLOSUM-62 matrix (18) gives a positive score are indicated by a ‘+’ symbol.

I think gringer is correct to say this is because of the BLOSUM62 matrix and a ‘+’ symbol inidcating a positive score for a mismatch but here is a little more explanation.
BLASTP looks for (nearly) identical short protein matches pairwise, called seeds, and then extends the seeds in each direction one residue at a time using arbitrary gap penalties for indels and a matrix for substitutions. Based on some thresholds and more, the highest scoring alignments are kept and ultimately reported.
The default choice of matrix for BLASTP is BLOSUM62, which looks like this. Obviously matches have high positive scores, but also mismatches have varying positive or negative scores.
If you look up your changes:
H -> Q; ; score 0
E -> K; +; score 1
Note that this matrix is symmetric so the scores for H -> Q == Q -> H and E -> K == K -> E.
In terms of why E -> K has a positive score, well not all amino acid changes are viewed equally. BLOSUM matrices are not based on any biological theory but instead use the log-odds ratio, a likelihood approach, that estimates from real sequence alignments how frequently different amino acid changes occur. The more often E -> K substitutions are observed in these alignments, the higher the substitution will score in the derived matrix.
Hopefully these scores are intupretable based on the biochemical properties of the amino acids. For example, R -> K scores 2, a  high score for a mismatch, but both amino acids are basic and so share a similar structure, presumably making a change more likely to occur and/or more likely to maintain the protein function and so not be harmful.
For your case, I don't have a good biological explanation for why E -> K scores positive, another author notes this is surprising, but ultimately the reason it is given a positive score is just because regardless of why, E -> K does apparently frequently occur, at least in the alignment data used to derive this scoring matrix. Incidentally, BLOSUM matrices are pretty ancient and based on few alignments by modern standards, but nonetheless appear to work quite well.


text:
A Database for Plant Metabolome or Secondary Metabolites?
I am searching for a plant metabolome or secondary metabolite database for a college bioinformatics project. Up until now, I have found these two databases, as I am giving their links below;

http://scbt.sastra.edu/pmdb/
https://www.plantcyc.org/https://www.plantcyc.org/

and searched for different scientific articles and reviews;

https://www.frontiersin.org/articles/10.3389/fpls.2013.00073/full
https://www.ncbi.nlm.nih.gov/pubmed/27258266

But only one and updated database which could contain all known plants and all of their known metabolites would be perfect. Could you suggest me any database or any other source with a high coverage for my project?

Can you use the plant reactome db? http://plantreactome.gramene.org/

Do you really need to be all plants, and all metabolites? We don't have even for humans all the metabolites! Let alone know all plants, and from all have all the metabolites. Welcome to the site! Nice first question

Thank you for your suggestion @b.nota, I will check the db you shared :)

Thank you for your feedback @Llopis, I have meant all "known" plants and metabolites, could you recommend me any db with a coverage as high as possible?

PhytoMine is an InterMine interface to data from Phytozome, the Plant Comparative Genomics portal of the Department of Energy's Joint Genome Institute, which hosts 81 assembled and annotation genomes, from 78 viridiplantae species.


text:
visualisation of genome alignment
I was asked to check the synteny of some genes in two genome assemblies of the same species (PacBio & Illumina). I was given two scaffolds couples (so 4 scaffolds total) in the Illumina genome and asked if they were assembled together in the PacBio assembly.
So, I proceeded with a DNA-DNA alignment with minimap2 of the two Illumina scaffolds on the PacBio whole genome. But I am not sure on how to directly use the sam file for the question I was asked. So I tried to visualize this alignment using Mauve. It worked for the first couple, but it failed (blank page) for the second couple of scaffolds.
Surely, Mauve is not the best option out, so can you suggest an alternative? This is a chordata genome.

After a bit of digging, finally found exactly what I was looking for. I'll share here in case someone else was looking for the same thing.
I used a standard Nucleotide BLAST (blastn), coupled with MUMmer. Short tutorial.
Assuming that xn & yn are the couple of scaffolds I want to check if they are assembled together in the Pacbio assembly (a multi-fasta file):
1) first, I blasted the pair of Illumina scaffolds (as separate queries) to the PacBio genome (as subject):
blastn -query ../illumina_asm/scaffold_x1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_x1; blastn -query ../illumina_asm/scaffold_y1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_y1
The candidates PacBio scaffolds shall appear in the top blastn hits in both out_scaffold_x1 & out_scaffold_y1.
Let's assume we found one common scaffold only - corresponding to the fasta entry pacbio_scaffold in the PacBio genome.
2) extract the sequence of the target scaffold from the PacBio assembly, using samtools faidx:
samtools faidx ../pacbio_asm/genome.fa; samtools faidx ../pacbio_asm/genome.fa 'pacbio_scaffold' > pacbio_scaffold.fa
3) concatenate the Illumina scaffolds to analyze:
cat ../illumina_asm/scaffold_x1.fa ../illumina_asm/scaffold_y1.fa > illumina_scaffolds.fa
4) compare the Illumina scaffolds to the PacBio assembly using nucmer from MUMmer:
nucmer --prefix=pacbio_scaffold pacbio_scaffold.fa illumina_scaffolds.fa
this will generate a pacbio_scaffold.delta file, used in the next step
5) filter the delta file to keep only the most meaningful alignments file using delta-filter from MUMmer (avoiding confusing messy plots later):
delta-filter -m pacbio_scaffold.delta > pacbio_scaffold.delta.m
6) use mummerplot from MUMmer to achieve a graph:
mummerplot --png pacbio_scaffold.delta.m -R pacbio_scaffold.fa -Q illumina_scaffolds.fa --prefix=pacbio_scaffold -large -layout
On the x-axis, the PacBio scaffold they have in common. On the y-axis, the two Illumina scaffolds (separated by the horizontal grey abline). Colors correspond to the orientation.
Example of Illumina scaffolds linked in the PacBio assembly:

If you prefer dotter plot: https://peerj.com/preprints/26567/. There is also minimap2 -DP ref.fa query.fa|miniasm/minidot - > dot.eps, though the visual is not as good as mummerplot.

Also: (1) MashMap - repo has a dotplot script: https://github.com/marbl/MashMap  & (2) Assemblytics - requiring a nucmer delta file: http://assemblytics.com

I think you could use circos for this. The configuration is kind of tedious but it gives nice results. You would basically need to retrieve alignments coordinates and parse them into this format.
scaf1 19720 31917 scaf2 28307 34227
scaf1 28307 34227 scaf2 19720 31917

Example figure:

nice thanks! I has also been advised to try MUMmer (http://mummer.sourceforge.net), for easy dotplots

I remember a few options for pairwise genome visualisation (which may be outdated BTW):

ACT
YASS
aliTV
LASTZ

Thanks Léo! Maybe also worth of notice: (1) Geopard for GUI-dotplot: http://cube.univie.ac.at/gepard; (2) Ribbon - and actually (3) Mauve could also be run from command-line: https://sourceforge.net/p/ngopt/wiki/How_To_Score_Genome_Assemblies_with_Mauve/
Enough material for next bioinfo-DEE meeting I guess ;)


text:
Complete Genomics Junction file for SV Detection
At present, I have vcf file for each patient generated from tumor data. I am using virtual normal correction method for somatic mutation detection because I don't have normal data. For SNP or small variants detection, I have already converted vcf into tsv format but for SV, I need complete genomics junction file or high confidence junction beta file. Basically this file is used to represent SVs in the subject sequence. Look at the following pics for better idea :

and the header of the file is
ID    LeftChr    LeftPosition    LeftStrand    LeftLength    RightChr    RightPosition    RightStrand    RightLength StrandConsistent    Interchromosomal    Distance    DiscordantMatePairAlignment    ...

This file is specially generated and used by the company "Complete Genomics". This file format is basically used in SV representation. I have no idea about the meaning of junction. For more information, look at Data type [http://www.bx.psu.edu/~giardine/tests/tmp/DataFileFormats112.pdf].
My question is how can I generate cg junction file from the given vcf and bam file.Can anybody tell me how to get it?
Thanks.
p.s. I'm aware that Complete Genomics was acquired by BGI in $117.6M. Now the paper I am trying to understand can be obtained from this github link because I don't have matched normal data. In order to identify SVs, I need cg junction file. I don't know how to get it. That's why I am looking here for some help. Thanks.


text:
Seurat installation fails
When I am trying to install Seurat in R markdown file:
install.packages("Seurat")

it fails with the following errors:
Installing package into ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3’
(as ‘lib’ is unspecified)
Warning in install.packages :
  dependency ‘proxy’ is not available
also installing the dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’

trying URL 'https://cran.rstudio.com/src/contrib/igraph_1.2.1.tar.gz'

Content type 'unknown' length 2703426 bytes (2.6 MB)

downloaded 2.6 MB

trying URL 'https://cran.rstudio.com/src/contrib/dtw_1.18-1.tar.gz'

Content type 'unknown' length 849256 bytes (829 KB) downloaded 829 KB trying URL 'https://cran.rstudio.com/src/contrib/diffusionMap_1.1-0.tar.gz'

Content type 'unknown' length 39130 bytes (38 KB)
downloaded 38 KB

trying URL 'https://cran.rstudio.com/src/contrib/Seurat_2.3.0.tar.gz'

Content type 'unknown' length 846065 bytes (826 KB)

downloaded 826 KB
* installing *source* package ‘igraph’ ...
** package ‘igraph’ successfully unpacked and MD5 sums checked
installing to /sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs
** R
** demo
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded

Error in dyn.load(file, DLLpath = DLLpath, ...) :    unable to load shared object '/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs/igraph.so': libicui18n.so.58: cannot open shared object file: No such file or directory

Error: loading failed Execution halted ERROR: loading failed
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph’ Warning in install.packages :   installation of package ‘igraph’ had non-zero exit status

ERROR: dependency ‘proxy’ is not available for package ‘dtw’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/dtw’ Warning in install.packages :   installation of package ‘dtw’ had non-zero exit status

ERROR: dependency ‘igraph’ is not available for package ‘diffusionMap’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/diffusionMap’ Warning in install.packages :   installation of package ‘diffusionMap’ had non-zero exit status

ERROR: dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’ are not available for package ‘Seurat’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/Seurat’ Warning in install.packages :   installation of package ‘Seurat’ had non-zero exit status

How could it be fixed? What to try?
Operating system is CentOS Linux release 7.3.1611 (Core)

Update

I also tried development version of Seurat but it failed too:
http://satijalab.org/seurat/devinstall.html

Update

I tried installing Seurat using conda install -c bioconda r-seurat and it is not working. It is giving the error:

-bash-4.2$ conda install -c bioconda r-seurat
  Fetching package metadata .............
  Solving package specifications:
PackageNotFoundError: Packages missing in current channels:

r-seurat -> r-cowplot
r-seurat -> r-diffusionmap
r-seurat -> r-dtw
r-seurat -> r-ggridges
r-seurat -> r-ica
r-seurat -> r-pbapply
r-seurat -> r-ranger
r-seurat -> r-rcppprogress
r-seurat -> r-sdmtools
r-seurat -> r-tclust

We have searched for the packages in the following channels:

https://conda.anaconda.org/bioconda/linux-64
https://conda.anaconda.org/bioconda/noarch
https://repo.continuum.io/pkgs/main/linux-64
https://repo.continuum.io/pkgs/main/noarch
https://repo.continuum.io/pkgs/free/linux-64
https://repo.continuum.io/pkgs/free/noarch
https://repo.continuum.io/pkgs/r/linux-64
https://repo.continuum.io/pkgs/r/noarch
https://repo.continuum.io/pkgs/pro/linux-64
https://repo.continuum.io/pkgs/pro/noarch

You seem to be missing the libicui18n.so library. Please [edit] your question and tell us what operating system you are using. If Linux, tell us what distribution.

I just updated. It is CentOS. I am working with rstudio on the cluster.

Ah, thanks. So what is the output of yum search libicui18n (if that returns nothing, try yum search libicui or even yum search icui until you get a result). Then, install the package. If that works, let me know and I can post answer (or just post it yourself, up to you).

The output in all the three cases is No matches found. Also, I do not have root privileges obviously, so maybe I should download this libicui18n.so library and install it locally somehow?

Hmm, I haven't used an rpm-based system in ages. You might need to search for libicui* or something. I did check CentOS's repository, and there is such a package. I don't remember how to install rpms locally though and I'm on mobile, so I can't help much. You could try asking how to find and install the relevant package on [unix.se].

As @terdon mentioned in the comments, you need to sudo yum install libicu on the relevant nodes (or have your admin do that for you). Alternatively, seurat can be installed via conda, which means you don't need root access.

Actually I tried installing Seurat using conda install -c bioconda r-seurat  and it is not working. It is giving PackageNotFoundError.

You need -c conda-forge too. You might then need to fix readline with conda upgrade -c conda-forge readline.

If you don't wish to use conda as your r package manager (because, for a start, it will install a new copy of R with its own library), you can install just ICU with conda conda install -c conda-forge icu, and then install seurat in R


text:
Protein expression units?
I'm starting for protein expression analysis.
In gene expression, we have TPM, RPKM etc. However, all I can see for protein expression is: "High", "Medium", "Low" and "Not detected".

https://www.proteinatlas.org/ENSG00000140326-CDAN1/tissue

Q: Is there a unit for measuring protein expression like gene expression? I've seen people talking about PE (protein existence). Can we use it for protein expression?

I would guess that the real unit for measuring protein expression would be pM (pico-moles) or a count of protein molecules.
I'm not sure we have any methods to calculate such quantities though.

The unit of measurement in protein expression is more or less "abundance". In mass spec this is derived from the area under the curve in the spectrum. In things like an ELISA this is derived from a color intensity.
Protein existence is a vague summarization of multiple data sources to try to quantify the likelihood of a protein existing. You can't use that for expression, though expression would influence the likelihood of existence.


text:
groupGO parameters explanation (ont and level)
I need to use groupGO in clusterProfiler to find functional profile of a list of genes, and I am having trouble finding out what some parameters of the function mean and which I should select for my specific case.
The function is defined as:
groupGO(gene, OrgDb, keyType = "ENTREZID", ont = "CC", level = 2,
  readable = FALSE)

I have mouse genes dataset. What I should select for ont and level parameters? What does ont types: MF, BP, CC mean? What is level? Where I could find this info? It is definitely not in the vignette...

That's funny, I asked the same question today on bioconductor support about the level argument. I can tell you what ont stands for, though. It is the sub-ontology, MF = molecular function, BP = biological process, and CC = cellular component. Read more about it at the GO consortium.

As mentioned in the comments MF is molecular function, BP is biological process and CC is cellular compartment. These are the 3 domains of the ontologies. The level refers to the level in the ontology graph. In the example in that link, pigmentation would be level 1, pigmentation during development level 2, regulation of pigmentation during development is level 3 and so on. You can find more discussion of this in this post on biostars from a year ago.
Regarding the exact settings you should use, I don't personally find CC to be an informative domain. The first couple levels are typically too generic to be useful and they'll include too many genes anyway for you to likely pick up any changes. I'd start around level 3 and see how things go (in other programs you would get all of the levels and domains at once).


text:
clusterProfiler, groupGO: meaning of GeneRatio
I am running groupGO function from clusterProfiler:
ggo <- groupGO(gene     = head(gene, 200),
           OrgDb    = org.Mm.eg.db,
           keyType = "ENSEMBL",
           ont      = "BP",
           level    = 2,
           readable = TRUE)

And getting the following output:
> head(ggo)
                   ID            Description Count GeneRatio
GO:0000003 GO:0000003           reproduction     5     5/200
GO:0008152 GO:0008152      metabolic process    40    40/200
GO:0001906 GO:0001906           cell killing     0     0/200
GO:0002376 GO:0002376  immune system process    20    20/200
GO:0006791 GO:0006791     sulfur utilization     0     0/200
GO:0006794 GO:0006794 phosphorus utilization     0     0/200

As can be seen cell killing, sulfur utilization and phosphorus utilization have zero genes. Then, why are they here at all? Should not the groupGO show only biological processes the fed in genes are related to? This seems very weird to me. What is the meaning of the GeneRatio column?

If you don't want terms with zero genes, just use ggo[ggo$Count>0]

As you seem to be learning, the clusterProfiler package is a bit "special". It doesn't filter the results you're seeing for relevance, it returns everything. The GeneRatio is the number of genes in your input list associated with the given GO term (e.g., 5) and then a / and the total number of input genes. This is generally not that useful. Since many things aren't documented, you may find the source code helpful.