text:
No variant found using GATK 4.0 HaplotypeCaller
I am doing variant calling on RNA-seq datasets from wheat which is hexaploid,the binary alignment (BAM) files were created using STAR version 2.6.0c and variant calling was done using GATK 4.0 HaplotypeCaller.The whole pipeline is as follows:
Alignment through STAR 2.6.0c :
STAR --runMode alignReads --runThreadN 20  --outSAMtype BAM  SortedByCoordinate --sjdbOverhang 100  --genomeDir wheat --readFilesIn G1_cleaned_R1.fastq G1_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3

Then AddOrReplaceReadGroups from picard was used with default settings:
 AddOrReplaceReadGroups \
   I=input.bam \
   O=Read_groups_added.bam \
   RGID=4 \
   RGLB=lib1 \
   RGPL=illumina \
   RGPU=unit1 \
   RGSM=20

Finally GATK 4.0 haplotypeCaller was applied:
gatk HaplotypeCaller --native-pair-hmm-threads 28  \
    -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf

The command ran successfully as shown below:
20:33:01.018 INFO  ProgressMeter - TGACv1_scaffold_470693_5DS:301            102.2              45108670         441320.0
20:33:11.039 INFO  ProgressMeter - TGACv1_scaffold_724195_U:301            102.4              45111270         440625.5
20:33:21.072 INFO  ProgressMeter - TGACv1_scaffold_724579_U:301            102.5              45113880         439932.4
20:33:31.108 INFO  ProgressMeter - TGACv1_scaffold_248869_3B:301            102.7              45116440         439240.9
20:33:41.143 INFO  ProgressMeter - TGACv1_scaffold_248918_3B:301            102.9              45118990         438551.6
20:33:51.150 INFO  ProgressMeter - TGACv1_scaffold_306137_4AL:301            103.0              45121550         437866.7
20:34:01.176 INFO  ProgressMeter - TGACv1_scaffold_193559_2DS:301            103.2              45124100         437182.5
20:34:03.134 INFO  HaplotypeCaller - 100112154 read(s) filtered by: ((((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter) AND WellformedReadFilter)
100112154 read(s) filtered by: (((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter)
100112154 read(s) filtered by: ((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter)
100112154 read(s) filtered by: (((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter)
100112154 read(s) filtered by: ((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter)
100112154 read(s) filtered by: (((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter)
100112154 read(s) filtered by: ((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter)
100112154 read(s) filtered by: (MappingQualityReadFilter AND MappingQualityAvailableReadFilter)
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

20:34:03.134 INFO  ProgressMeter - TGACv1_scaffold_725922_U:301            103.2              45124592         437049.1
20:34:03.134 INFO  ProgressMeter - Traversal complete. Processed 45124592 total regions in 103.2 minutes.
20:34:04.749 INFO  VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
20:34:04.749 INFO  PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.0
20:34:04.749 INFO  SmithWatermanAligner - Total compute time in java Smith-Waterman : 0.00 sec
20:34:04.749 INFO  HaplotypeCaller - Shutting down engine
[July 3, 2018 8:34:04 PM PKT]   org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 103.99 minutes.
Runtime.totalMemory()=4789370880

but still the generated vcf file i.e final.vcf does not contain any variants as shown below:
[dellemc@localhost bilal]$ tail final.vcf
##contig=<ID=TGACv1_scaffold_725916_U,length=500>
##contig=<ID=TGACv1_scaffold_725917_U,length=500>
##contig=<ID=TGACv1_scaffold_725918_U,length=500>
##contig=<ID=TGACv1_scaffold_725919_U,length=500>
##contig=<ID=TGACv1_scaffold_725920_U,length=500>
##contig=<ID=TGACv1_scaffold_725921_U,length=500>
##contig=<ID=TGACv1_scaffold_725922_U,length=500>
##contig=<ID=TGACv1_scaffold_731344_7BL,length=500>
##source=HaplotypeCaller
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  20

Is there any mistake in the haplotypeCaller command or any argument is missing ?

All those compute times of 0.0 seem a little suspicious.

I have edited the question to add more details.

@swbarnes could you tell about the type of suspicions?

There are steps of the 'best practices' workflow missing, which I think is leading to the majority of reads being filtered by the haplotype caller, notice these lines in the output:
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

That's >100m reads being filtered by the HC.
Specifically, you're missing SplitNCigarReads, which separates reads that bridge a splice junction (padded by Ns in the STAR alignment) and reassigns the mapping quality. At the moment this element of the workflow is not validated or available in GATK4 (RNA-Seq variant calling is largely abandonware in newer versions of GATK, AFAIK), but in the best practices WDL utilises GATK3.5 - see https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl, line 440-477.
There are other elements of the recommended workflow missing too, but I think this is the critical bit. I won't harp on here about the unsuitability of RNA-Seq data for variant calling, as I'm all too aware there often isn't much option.

What is meant by the option --limitOutSJcollapsed in the STAR command given in the link you suggested https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/a7a890b38af1309edc040ac89814df37a57f77b7/rna-germline-variant-calling.wdl#L345

A 'collapsed' junction is what gets written to the SJ.out.tab file during a run of STAR (the collapsing here refers to the removal of duplicate junctions). This option merely limits the maximum number of these junctions that are written. It defaults to 1,000,000.

As well as the issues @sjcockell mentions, you have another potential problem you've got a hexaploid organism, by default the HaplotypeCaller (HC) is hardwired to a ploidy of 2, you need to set this to 6 to match your organisms ploidy using -ploidy 6 as an argument.  Since the ploidy is a key term in the underlying Bayesian stats the HC employs (see this link for an explanation) it's likely that some of the lower allelic frequencies that hexaploidy would generate will mean quite a few variants won't be called as they'd be statistically unlikely for a diploid organism.

You might also want to consider the GATK support forum too, as they are very helpful on GATK matters: https://gatkforums.broadinstitute.org/gatk/categories/gatk-support-forum


text:
plotting two heatmaps with the same order of genes
I expected the same pattern but here I am not able to compare the patterns as the order of genes does not seem the same. I mean how I can have two heat maps on which the order of genes are the same so I woulds be able to compare the block of yellow and blue colors in two heat maps.
I have used this code
require("RColorBrewer")
myCol <- colorRampPalette(c("dodgerblue", "black", "yellow"))(100)
myBreaks <- seq(-2, 2, length.out=101)
heat <- t(scale(t(sc_DEGG)))
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete")
heatmap.2(
    heat,
    Rowv=as.dendrogram(hr),
    Colv=as.dendrogram(hc),
    col=myCol,
    breaks=myBreaks,
    main="Title",
    key=T,
    keysize=1.0,
    scale="none",
    density.info="none",
    reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
    trace="none",
    cexRow=0.2,
    cexCol=0.8,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

I  have made a pair of heat maps based on your code in complexheatmap package. The right heat map is lighter than the left one. I don't know which one has been clustered based on the another one. How I can make the right one darker so that both  seem more similar??

ComplexHeatmap is built for plotting side-by-side heat maps with the same clustering - you use the + notation, similar to ggplot2. To use the same example data as @b.nota:
library(ComplexHeatmap)
# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))
rownames(m) <- letters[1:26]

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))
rownames(m2) <- letters[1:26]

Heatmap(m) + Heatmap(m2)

And the results:

There are lots of examples in the extensive documentation.
EDIT:
Since you have a few genes with high values in the 'T' matrix, they are washing out the colour scheme for the rest of the heatmap. You could set the max to be lower, and accept that those genes will saturate - the resulting heatmap will be less 'honest' about the data, but might capture what you want to see. To change colour scheme requires a call to the colorRamp2 function of the circlize package.
library(circlize)
mycol <- colorRamp2(c(-2,0,2), c("dodgerblue", "black", "yellow"))
Heatmap(m, col=mycol, cluster_columns = FALSE) +
    Heatmap(m2, col=mycol, cluster_columns = FALSE)

This gives:

Great solution!

You expect that the order of genes is the same for both heatmaps, but the chance is small that the 4000 genes will exactly be clustered in the same order.
What you can do, is use the order of genes (rows) from one heatmap and create a new heatmap with the other data set. You'll lose the tree of course because it will not be clustering anymore.
Here a small example:
library(gplots)

# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))

rownames(m) <- letters[1:26]

hm <- heatmap.2(
    m,
    main="m",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))

rownames(m2) <- letters[1:26]

hm2 <- heatmap.2(
    m2,
    main="m2",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix with row order of first
reorder_m2 <- m2[rev(hm$rowInd),]

hm2_m <- heatmap.2(
    reorder_m2,
    main="m2 in m's order",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Rowv = FALSE,
    Colv = FALSE,
    dendrogram="none",
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

You have 2 options provided that both datasets have the same genes (rows):

Keep the original order of rows in the heatmap. If they’re in the same order in the dataset then they will be in the heatmap. You can suppress reordering with any of the following arguments to heatmap.2

Rowv = NULL
dendrogram = “col” or dendrogram = “none”
reorderfun = function(x) return(x)

You can compute the order of the genes (rows) in advance and pass them to the heatmap. The code shown in the question does this already with Rowv=as.dendrogram(hr). To have the same order in a second heatmap, you have to pass the same same dendrogram hr to both heatmaps. This can be computed from either dataset or another set of samples altogether.

Another option is to compute separate dendrograms (from hclust of each dataset) for the samples and combine them using merge. You can combine the datasets with cbind and plot the datasets together as one heatmap (with 2 sets of samples clustered separately) by passing the merged dendrogram to Colv:
#compute same dendrograms separately
hc1 <- hclust(as.dist(t(data1), method="euclidean")), method="complete")
hr2 <- hclust(as.dist(t(data2), method="euclidean")), method="complete")
mergedcols <- merge(as.dendrogram(hr1), as.dendrogram(hr2))
mergeddata <- cbind(data1, data2)
#compute gene dendrogram across all samples
hr <- hclust(as.dist(1-cor(t(mergeddata), method="pearson")), method="complete")
heatmap.2(
mergeddata,
Rowv=as.dendrogram(hr),
Colv=as.dendrogram(mergedcols),
col=myCol,
breaks=myBreaks,
main="Title",
key=T,
keysize=1.0,
scale="none",
density.info="none",
reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
trace="none",
cexRow=0.2,
cexCol=0.8,
distfun=function(x) dist(x, method="euclidean"),
hclustfun=function(x) hclust(x, method="ward.D2")
)

You should be able to show the division between these samples with the colsep argument.

source: Stack Exchange

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

sort: votes


text:
Meaning of category in MutSigCV
I am trying to understand the working of MutSigCV. But I am not sure about some terms used in MAF file format which is category.
According to MutSigCV, every mutation can be divided into the following categories:

CpG transitions
CpG transversions
C: G transitions
C: G transversions
A: T transitions
A: T transversions
null+indel mutations

I am very much clear about some of this categories like C: G transitions and A: T transitions but for rest, I am unable to make any solid conceptual background. Can anyone suggest me the proper definition of these categories of mutations with proper examples?
Thanks.

Definition
MutSig stands for "Mutation Significance" and is defined here

Did you try asking on the help forum? Or your problem is what is each mutation?

@Llopis: Not yet. I even didn't know about that. I know only GATK forum and this website. I have googled these terms extensively but didn't get even proper definitions. That's why I started the bounty on this question.

CpG transversions are the same as C: G transversions, except that they are followed by a G which means they are potential targets for DNA methylation. Methylcytosine can undergo deamination to Thymine, which makes them mutation hotspots

indel = insertion/deletion. I have no idea what null means in this context

@heathobrien Those comments are an answer! With a link to a source you could get the bounty

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


text:
Find paralogs in a draft genome
We generated a (diploid, chordata, highly heterozgous) genome using PacBio and we wanted to see whether it contains lineage-specific duplications (paralogs, basically). The genome is not in Ensembl yet.
The only data we have at the moment are:

genome
transcript annotation
RNAseq

We found some methods from papers:

use Blast
detect segmental duplications in complete genomes with SDDdetector
detecting putative recent segmental duplications or diploid homozygous organisms based on NGS data DuplicationDetector
(I just came up with this, but should work) map back the reads to the assembly & analyze the read depth to detect duplicated portions

I'll gladly take advice.

Chordata...highly heterozygous... Do I sense a tunicate?

nope, a lancelet ;)

If there are other genomes or transcriptomes close to your taxon, you might try to cluster homologs and infer gene trees from homologs groups to detect duplications.

thanks @NatWH. yes I have a genome/transcriptome (the stats are less good than our assembly, but should not be a terrible issue)

This is really tough to do with highly heterozygous animals. What are your genome assembly stats? Specifically, what is your number of contigs, scaffolds, assembly size, and the N50?
If you have proximity ligation data it will be easier to determine if potential paralogs are truly from different regions of the genome or are just from homologous locations on sister chromatids that ended up in uncollapsed in your final assembly.
If I were you and had a good genome, I would start by annotating gene models using RNAseq data, then blasting each gene in the model against all genes in the model. This would give you potential paralogs to start looking at more closely.

Thanks. About my PacBio genome: haplo size ~550Mbp (canu + purge haplotigs),  ~1000 contigs, ~1.6 Mbp N50, not yet scaffolded but not sure it's relevant. No proximity ligation data, unfortunately - but I was thinking of checking reliability of paralogs by making sure they are found in an alternative assembly we have (same species, but Illumina, so quite fragmented). I'll test what you suggest in 3rd paragraph too, thanks for the input!

Those are great numbers - I'd make sure that you polish with pilon then racon to remove indels before mapping RNAseq reads and generating gene models.

also @conchoecia  can you suggest an alternative to pilon? it is apparently "not currently tuned to the error model of raw PacBio reads, and their use may introduce false corrections" (https://github.com/broadinstitute/pilon/wiki/Methods-of-Operation)

Ah, for pilon I mean use Illumina data. You might use Arrow (pb reads) then pilon (Illumina reads).

oh ok! I don't have Illumina data yet. I'll try with Arrow (I tried with Racon already but only seemed to make the assembly worse -.- not sure why).


text:
Variability in genes across platforms
I have a matrix of gene expression for single cell (cells in column and genes in row) in 9 time points with about to 200 cells for each time point. I have also a matrix of gene expression in bulk RNA-seq with the same time point and 2 replications.
I have pooled reads from cells in each time point to have a simulated bulk RNA-seq. now I have two matrices; Real bulk RNA-seq in 9 time points and two replication and another matrix a simulated bulk RNA-seq with 9 time points and two replications. I have plotted a heat map of them to show the pattern of gene expression is similar in both matrices so single cell RNA-seq is good enough. Actually I have used the real bulk RNA-seq in a similar study as a gold standard for checking my single cell RNA-seq. My problem is this the pattern of gene expression in 0 time point is very different in two datasets. I should note that 0 time point comes from Fluidigm and another 8 time points come from icell8. I should be able TO adjust this heat map to look similar enough.
So I thought to show similarity by plotting cumulative fraction of genes detectable in varying fractions of cells (lets say 20%, 40%, ...), segregated by bulk RNA-seq expression quartile. Something like a plot that x axis is percent of cells and y axis is percent of genes but I think I should solve the fundamental problem in a heat map.
I have single cell (matrix with T) and bulk RNA-seq (matrix with R). Whatever I am plotting, time point 0 in single cell RNA-seq seems more yellow and dose not allow a good similarity between single cell and bulk RNA-seq. How I can make zero time point less yellow so another time points become more yellow?

How do you rank each cell by number of genes expressed?

Thank you, actually I am not sure, if I have FPKM values for bulk RNA-seq, I don't   know how to show how many cells express a quartile of genes in bulk RNA-seq

Sorry, I think I misunderstood the question: what has to do the bulk RNA-seq with identifying the variability of genes in individual cells? What do you try to show ? You expect that this graph (however it is made) will help you to find/learn that ...?

You right the question is very ambiguous,  I want to check similarity of two sets of samples

First make sure what you want help with in this question. Then please [edit] the question with more details. Do you want to compare the bulk RNA to the scRNA? Or you want to check the similarity at certain timepoints or between replications?

I'm kind of surprised that your method of comparing data from two totally different sample prep methods worked at all, not that it differs significantly at only one timepoint.

What do you want? A method to normalize between technologies so that your time point 0 looks more similar between scRNA and bulk or plot the similarity  between bulk and scRNA? Also I assume that this data wasn't corrected by batch/technology effect or it is?

I have not done batch correction. I only have pooled 200 cells in each time point to simulate a bulk RNA-seq and I hoped when comparing that with real RNA-seq from the same experiment the data look similar. But likely 0 time point is very different coming from two sequencing methods. after pooling my cells I normalised them with deseq2 like I have done for real bulk RNA-seq. for 8 time points looks good but for zero as this platform gives more read counts. how I manipulate my matrices to make a more similar heat maps or how could I remove difference between 0 time point and another time points

Sorry you are right about different sample prep but when I am calculation Pearson correlation between each time point I am seeing a high correlation between data. If this high correlation is by chance? correlation after log transforming both of datasets


text:
pysam pileup: what reads appear in the pileup?
Cross-posted on Biostars (with no answers currently). I hope that's OK.

I thought that iff a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.
This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.
What's wrong in my understanding?

Updated code (and updated output) moved to Gist.

Edit: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.
E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-<  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69

I guess this explains why the read appears in more than 150 pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in 198 = (48 + 6 + 144) pileup columns, at each position where one of them maps?

Update: I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case pileup() was filtering reads with low base quality, I also calculated how many positions the read has baseQ >= 13. It's 150 positions, but this read still appears in pileup columns for only 121 positions.

Can you include the raw sam lines for the read so we can determine if the issue is in pysam or your code.  That would give us truth.

@Bioathlete please see above

The two entries in the sam file represent mate pairs for the given ID.  You can tell the difference based on the sam flags.  Picard has a nice tool to determine the meaning of the flags.
You are missing coverage within the first read.  I think that it is related to the low quality bases.  samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.

Thanks for the tip about distinguishing a read and its mate by QNAME and flag—that's helpful. The read with flag 83, though, has baseQ > 13 at 150 positions, yet appears in the pileup columns for only 121 positions.

Ah but BAQ ("Base Alignment Quality"), not just regular base quality, also has to be above 13. Comparing samtools mpileup and samtools mpileup -Q0 shows that that's why reads are being filtered out.


text:
tRNAscan-SE error: FATAL: Unable to find /usr/local/bin/cmsearch executable
I have downloaded tRNAscan-SE from here. After decompressing and untaring the file, I installed it using:
./configure
make
make install

When I type tRNAscan-SE --help, I get the help page:
tRNAscan-SE 2.0 (December 2017)

FATAL: No sequence file(s) specified.

Usage: tRNAscan-SE [-options] <FASTA file(s)>

  Scan a sequence file for tRNAs
          -- default: use Infernal & tRNA covariance models
             with eukaryotic sequences
             (use -B, -A, -M, -O or -G to scan other types of sequences)

Basic Options
  -E         : search for eukaryotic tRNAs (default)
  -B         : search for bacterial tRNAs
  -A         : search for archaeal tRNAs
  -M <model> : search for mitochondrial tRNAs
                 options: mammal, vert
  -O         : search for other organellar tRNAs
  -G         : use general tRNA model (cytoslic tRNAs from all 3 domains included)
  -L         : search using the legacy method (tRNAscan, EufindtRNA, and COVE)
                 use with -E, -B, -A, -O, or -G
  -I         : search using Infernal (default)
                 use with -E, -B, -A, -O, or -G
  -o <file>  : save final results in <file>
  -f <file>  : save tRNA secondary structures to <file>
  -m <file>  : save statistics summary for run in <file>
               (speed, # tRNAs found in each part of search, etc)
  -H         : show both primary and secondary structure components to
               covariance model bit scores
  -q         : quiet mode (credits & run option selections suppressed)

  -h         : print full list (long) of available options

I then tried to use it with tRNAscan-SE c_elegans.PRJNA13758.WS263.genomic.fa but it gave the following error:
FATAL: Unable to find /usr/local/bin/cmsearch executable
After googling cmsearch, it seems that is is part of the Infernal package. So, I installed infernal with conda using conda install infernal.
I pulled the help from cmsearch successfully:
cmsearch -h
# cmsearch :: search CM(s) against a sequence database
# INFERNAL 1.1.2 (July 2016)
# Copyright (C) 2016 Howard Hughes Medical Institute.
# Freely distributed under a BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: cmsearch [options] <cmfile> <seqdb>

Basic options:
  -h        : show brief help on version and usage
  -g        : configure CM for glocal alignment [default: local]
  -Z <x>    : set search space size in *Mb* to <x> for E-value calculations  (x>0)
  --devhelp : show list of otherwise hidden developer/expert options

Options directing output:
  -o <f>       : direct output to file <f>, not stdout
  -A <f>       : save multiple alignment of all significant hits to file <s>
  --tblout <f> : save parseable table of hits to file <s>
  --acc        : prefer accessions over names in output
  --noali      : don't output alignments, so output is smaller
  --notextw    : unlimit ASCII text output line width
  --textw <n>  : set max width of ASCII text output lines  [120]  (n>=120)
  --verbose    : report extra information; mainly useful for debugging

Options controlling reporting thresholds:
  -E <x> : report sequences <= this E-value threshold in output  [10.0]  (x>0)
  -T <x> : report sequences >= this score threshold in output

Options controlling inclusion (significance) thresholds:
  --incE <x> : consider sequences <= this E-value threshold as significant  [0.01]
  --incT <x> : consider sequences >= this score threshold as significant

Options controlling model-specific reporting thresholds:
  --cut_ga : use CM's GA gathering cutoffs as reporting thresholds
  --cut_nc : use CM's NC noise cutoffs as reporting thresholds
  --cut_tc : use CM's TC trusted cutoffs as reporting thresholds

Options controlling acceleration heuristics*:
  --max      : turn all heuristic filters off (slow)
  --nohmm    : skip all HMM filter stages, use only CM (slow)
  --mid      : skip first two HMM filter stages (SSV & Vit)
  --default  : default: run search space size-dependent pipeline  [default]
  --rfam     : set heuristic filters at Rfam-level (fast)
  --hmmonly  : use HMM only, don't use a CM at all
  --FZ <x>   : set filters to defaults used for a search space of size <x> Mb
  --Fmid <x> : with --mid, set P-value threshold for HMM stages to <x>  [0.02]

Other options*:
  --notrunc     : do not allow truncated hits at sequence termini
  --anytrunc    : allow full and truncated hits anywhere within sequences
  --nonull3     : turn off the NULL3 post hoc additional null model
  --mxsize <x>  : set max allowed alnment mx size to <x> Mb [df: autodetermined]
  --smxsize <x> : set max allowed size of search DP matrices to <x> Mb  [128.]
  --cyk         : use scanning CM CYK algorithm, not Inside in final stage
  --acyk        : align hits with CYK, not optimal accuracy
  --wcx <x>     : set W (expected max hit len) as <x> * cm->clen (model len)
  --toponly     : only search the top strand
  --bottomonly  : only search the bottom strand
  --tformat <s> : assert target <seqdb> is in format <s>: no autodetection
  --cpu <n>     : number of parallel CPU workers to use for multithreads

*Use --devhelp to show additional expert options.

I tried the tRNAscan-SE command again but it failed with the same error message.
So, tRNAscan-SE is looking for cmsearch in the wrong place. How can I make tRNAscan-SE use cmsearch successfully?
Cheers.
Update
Forgot to tell you I am working on Mac OS X and that cmsearch is installed here:
/Users/user/anaconda/envs/python3env/bin/cmsearch
and my python3env anaconda environment is active.

Please [edit] your question and tell us i) what operating system you are using and ii) assuming it is some sort of *nix, the output of which cmsearch

An easy way to build tRNAscan and their dependencies "all Infernal packages" without confusion is to run
conda install -c bioconda trnascan-se

You may try to set the correct path through the infernal_dir variable in tRNAscan-SE.conf.src. You might need to do this before the installation process.
Original file starts with:
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: {bin_dir}

Try the following:
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: "/Users/user/anaconda/envs/python3env/bin"

I've done what you said and it still works but how do I know it's not because my symbolic links are still active? Do you know how I can erase them?

@charlesdarwin You can remove them by  (carefully) performingsudo rm /usr/local/bin/cmsearch and sudo rm /usr/local/bin/cmscan.

tRNAscan-SE 2.0 requires separate installation of Infernal 1.1.2.
The source code and prebuilt binaries can be obtained at
http://eddylab.org/infernal/.
The installation directory of Infernal should be the same as the one
where tRNAscan-SE 2.0 is installed.
When install infernal, run the configure command with the option --prefix <your tRNAscan-SE-2.0 install path>.
If you didn't add option --prefix when installing tRNAscan-SE-2.0, then you should install infernal without option --prefix.

A much cleaner solution is given in @bli's answer.

The error message is telling you that tRNAScan-SE is looking for the cmsearch executable at /usr/local/bin/cmsearch. However, you installed it though conda which apparently placed it at /Users/user/anaconda/envs/python3env/bin/cmsearch.
A simple, if not very elegant solution, would be to make a link to it:
sudo ln -s /Users/user/anaconda/envs/python3env/bin/cmsearch /usr/local/bin/

had to do the same with cmscan and now the code is running.

what does this soft link do?

@charlesdarwin This enables the program to access the program from the place it looks at by default. There might be a cleaner way to do: Changing the default (see my untested answer).

@charlesdarwin as bli said, that just creates a "fake file" at /usr/local/bin/ which points to /Users/user/anaconda/envs/python3env/bin/cmsearch, so that when something tries to open /usr/local/bin/cmsearch, it will actually run /Users/user/anaconda/envs/python3env/bin/cmsearch. But bli's approach of changing the configuration file is a much cleaner solution. Please use that instead.


text:
finding RNA-protein physical interaction
We all know that for physical protein-protein interaction, we need to find the distance between residues from PDB file of that interaction (finding distance between carbon alpha, carbon beta or centroid of two residues in PDB data of two proteins). I need to investigate the physical interaction between a protein and an RNA and I have PDB ID for this interaction but I don't exactly know to find the distance between which parts of residues (calpha, cbeta, centroid and etc.) with which parts of nucleotides (bases, phosphate group, or ribose backbone) from a PDB ID. Is there anyone who knows what can I do for this purpose?

I have reservations about the "We all know..." part. I'd use distance between the interacting atoms rather than between arbitrary points of the residue.

The question seems to be a follow-up of https://bioinformatics.stackexchange.com/questions/860/best-distance-parameter-for-estimating-physical-interaction-between-residues-in

@marcin, that question was about residue-residue interactions but this question is about nucleotide-residue interactions. For my special purposes I didn't need to find the distance between all atoms but I just needed to find the distance between arbitarary points of the residues and here in this question between residue and nucleotides

If you were to treat the nucleotide like a ligand you could follow these steps on pymol, and find the measurements you desire.
Calculating distances to centroids is a little bit more complex, as you have to make a pseudo atom to measure against.


text:
gene-level versus transcript-level analysis
Traditionally, RNA-seq data was quantified on gene level. Newer methods quantify on transcript/isoform level. For example, Kallisto only outputs transcript-level abundances. From the DESeq2 vignette:

A newer and recommended pipeline is to use fast transcript abundance
  quantifiers upstream of DESeq2, and then to create gene-level count
  matrices for use with DESeq2 by importing the quantification data
  using the tximport package.

Is it just for consistency and/or simplicity that the values are converted to gene level? Would it better to proceed with transcript-level table or does that violate some DESeq2 (or similar tools) assumptions?

Comments on this answer may help a little: https://bioinformatics.stackexchange.com/a/4508/2841

There are a variety of reasons people use gene-level quantitations.

Transcript-level differences are difficult to biologically interpret. Let's be honest, few groups are likely to put in the work required to determine what these might mean. Most genes have at least some level of characterization, so it's much more biologically tractable to think in terms of gene-level changes. Granted, "genes" are largely artificial constructs and in an ideal world one would always think in terms of transcripts, since they're what's really changing. But we don't live in an ideal world.
Statistical power increases with counts, so if you're splitting your counts among a handful of isoforms you could be missing something. In practice this isn't a big problem, since there are typically obvious major isoforms, but that won't always be the case.
"Mo transcripts mo problems", to channel the bioinformatics notorious B.I.G. It's common practice for people to use a subset of transcripts when using tools like Kallisto and Salmon, since the results usually end up being better. But which subset of transcripts should you use? Those in refseq? Gencode basic? A subset of one of those? Do we really understand all the consequences of choosing one over the others? (No, we don't) We all have our own preferences, but this isn't a solved problem.

(4) it is not at all clear that using transcripts rather than genes would be more informative in the vast majority of cases, even if we had perfect data (there are some people who disagree but most people seem to agree).

To add to the list that Devon Ryan outlined (or perhaps to elaborate on point 2?):
Although Salmon/Kallisto/RSEM are the more accurate in their transcript quantification than the methods they superseded, transcript level quantification is still not as accurate as gene level quantification (see the tximport paper, which would also be the tool i'd recommend for getting gene level estimates from the transcript-level output of these tools).
That said, the tximport authors do claim that considering the transcript-level quantification in your gene-level analysis improves the accuracy of the gene-level analysis.


text:
Associating SNP and GENE
Assuming I have SNPs data using hg19, how can I know which SNP belongs into which Gene?
The data looks like:
chr10_103577643
chr10_124712463

and so on.
I want to add a column of Gene, which would tell to which Gene the SNP belongs.
The file is not a VCF file but a Txt file.

Please [edit] your question and show us an example of your data. What do you have? A list of rsIDs? A VCF file? Some other format? How many SNPs?

Just annotate the vcf file.

Use the Variant Effect Predictor(VEP)

Please don't just post links as answers. It would be great if you could add a short explanation of what VEP is, how to use it etc. Otherwise, this is really a comment and not an answer and would become useless in the (admittedly unlikely) case where the VEP were to be discontinued or its URL changed. For instance, pasting chr10_103577643 into VEP doesn't work, so you could explain what format the data would have to be in.

You can use snpEff or snpSIFT java annotator also for hg19.

Please don't just post links as answers. It would be great if you could add a short explanation of these tools do, how to use them etc. Otherwise, this is really a comment and not an answer and would become useless if the tools were to be discontinued or their URLs changed.


text:
What is tRNAscan-SE bit score?
I have used the command-line version of tRNAscan-SE by Lowe et al. (1997) on the reference genome sequence of C. elegans to predict tRNAs. The programme outputs a table with tRNA locations, amino acid, anti codon, intro locations, tRNA type (functional or pseudogene) along with a score. The score in my dataset ranges from 20.0 to 107.4. After reading the article Lowe et al (2016) describing this tool, I still do not understand how this score comes about and what it means. I understand that:

the higher the score, the more likely the tRNA is assigned to the correct isotype (amino acid)
the lower the score, the more likely the tRNA is a pseudogene

My questions are:

How is the score calculated?
What does the score mean?

The authors of the tool state that is is 99% sensitive and detects 1 false tRNA in 15 GB of DNA. Does this mean that I do not need to filter my results? What about the confidence in the type of tRNA and whether it is functional or pseudogenised?
Lowe et al (1997) tRNAscan-SE:  a  program  for  improved  detection  of transfer  RNA  genes  in  genomic  sequence. Nuc. Acids 25
Lowe et al (2016) tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes. Nuc Acids Web Server issue: W54–W57.


text:
How many reads do I need to cover the entire genome?
Suppose my genome is 3 million bases and that my reads are 100 nucleotide long. I need to know how many reads I need to cover the entire genome.
I start from using the equation $C = \frac{N \cdot L}{G}$ where C is the coverage, N the number of reads, L the length of a read and G the length of the haploid genome. I also know that the average coverage follows a Poisson distribution, so I need this for C.
On my slides, in fact, I have the following:
$$
C = \frac{N \cdot L}{G} \approx ln \left( \frac{G}{L \cdot \epsilon}\right),
$$
if we assume that this is the average coverage needed to cover the entire genome with probability $1-\epsilon$.
The problem is that I do not understand how to arrive at the expression on the right. I thought that if we want to be sure to have a coverage of at least 1 on average, we would have had to compute $1 - P(X=0) = 1 - e^{-\frac{N \cdot L}{G}} = 0.99$, if we take $\epsilon=0.01$. But I do not get the same result as with the equation above.

Let's see what your slides are claiming $\epsilon$ to be, using:
$$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$
We can rearrange to get
$$ \epsilon \approx \frac{G}{L} \times e^{-C} $$
Now, what does it mean for an entire genome to be covered? Having a clear understanding of this definition is crucial, and your misunderstanding seems to partially stem from this.
I claim it's natural to consider a genome as consisting of $G/L$ bins. In this setting, we have $N$ balls that each randomly goes into a bin. Then, the entire genome is covered when all bins contain at least $1$ ball/read.
From the Poisson distribution, we know that P(a bin has 0 balls) = $e^{-C}$.
Then, we define $\epsilon$ to be the probability that at least one of the $G/L$ bins has 0 balls. Given the way we've modeled the problem, this is exactly equal to
$$ \epsilon = 1 - (1 - e^{-C})^{G/L}. $$
From above, your slides are claiming $ \epsilon \approx \frac{G}{L} \times e^{-C}$. This might be strange, as it appears to be equal to the expected number of bins with 0 balls, instead of the probability that there is at least 1 bin with 0 balls. Furthermore, it's strange because it contains only one exponential, compared to our exact derivation which has a double exponential.
The mismatch in the number of exponential is a hint that another mathematical approximation is at play. We can use a first-order approximation of the binomial series to obtain
$$ (1-e^{-C})^{G/L} \approx 1 - (G/L)e^{-C} $$ which holds particularly when $|e^{-C}| \times \frac{G}{L} << 1$.
In your case, $\frac{G}{L} = 30,000$, so the approximation only really starts getting accurate when $C > 15$ or so.
Following through with our approximation, we get
$$ \epsilon \approx 1 - (1 - (G/L) e^{-C}) $$
$$ \epsilon \approx \frac{G}{L} \times e^{-C}. $$


text:
Where to download baseline/average gene expression level of all human coding genes?
I am looking for the most appropriate dataset for downloading baseline gene expression level across all human coding genes during development. I am aware that EMBL Expression Atlas is one of the resources that provide such information, but I am ideally looking for species averaged (and probably scaled) expression data for all human coding genes during early stages of development!
Any insight is much appreciated!

There is not going to be a public database with this data. Apart from anything else, generating data for early human development is difficult and ethically tricky. Also, most people who care about development would probably think that a whole body average would be meaningless.
The best you could do would be to select tissues you were interested in from something like Expression Atlas and calculate the averages yourself.


text:
RNA seq fasta file annotation from alignment to reference matches
I've got a fasta file with some RNA seq data and another csv file with the output from plast where I've aligned it to a reference using plastn. I'm struggling with figuring out a command to append my original RNA seq data fasta headers with the annotations.
I've looked at awk and sed but can't wrap my head around how to use them in this instance.
CSV output example:
Contig11    ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA   78.95   95  20  0   371 465 4627    4721    1e-14   82.4    75  625 1   0   15.2    0   5692    1   0   15.2    0

And example fasta contig:
>Contig11
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

How would I be able to parse through the csv to find a matching contig with the fasta file, and append the second csv column to the fasta header. Ideally to produce:
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA....

Cheers for any help!

What format is your "csv" in? It isn't a "real" csv (that means "comma separated file" where columns are defined by commas), so what is the column separator? Is it a space or a tab?

An awk approach (assuming your input "csv" file is space-separated):
$ awk '(NR==FNR){a[">"$1]=$2; next;} (/>/){$1=$1" "a[$1]}1;' foo.csv foo.fa
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

Explanation

NR==FNR : FNR is the line number of the current input file and NR is the line number of all input. When processing more than one file, the two will only be equal while reading the first file.
(NR==FNR){a[">"$1]=$2; next;} : if we're reading the first file, the csv, save the second field as the value in the associative array a, whose indices are a > followed by the 1st field. I use the > only to simplify matching it later in the fasta file. The next moves on to the next line so the rest of the command isn't run for the 1st file.
(/>/){$1=$1" "a[$1]}: Now we're in the 2nd file. If the line starts with a >, then set the 1st field to be itself, a space and then whatever was stored in the array a for this 1st field.
1; : this is awk shorthand for "print this line".

You could expand the above into:
awk '{
        if(NR==FNR){
            a[">"$1]=$2;
        }
        else{
            if(/>/){
                $1=$1" "a[$1]
            }
        }
        print;
      }' foo.csv foo.fa

If your files are small enough to be sorted easily, you can use the FastaToTbl and TblToFasta scripts I have posted before and then join:
$ join -o 1.1,2.2,1.2 <(FastaToTbl foo.fa | sort ) <(sort  foo.csv ) | TblToFasta
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
>Contig12 mamamama|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
AKHSALKJHDLKAJBDLKAJDBADADKJALDKJBAD


text:
Selection of differential expressed genes
I'm working with RNA-seq data. I have 40 tumor samples and 5 Normal samples. Differential analysis with Deseq2 based on Fold change > 1.2 and alpha < 0.05 gave very low number of differentially expressed genes. Only 2 upregulated genes.
res <- results(dds, lfcThreshold = log2(1.2), alpha = 0.05)

To get more number of differential expressed genes I have few questions now

Instead of FDR < 0.05 can I use FDR < 0.1 (or) FDR < 0.5. Will there be any low confidence with this?
Can I select differential expressed genes only based on FDR < 0.05 without any fold change cutoff?
As I get very low number of DEGs with FDR < 0.05 & Foldchange > 1.2, Can I select DEGs based on Foldchange > 1.2 and p.value < 0.01 or 0.05 ?

PCA of the samples looks like this:

(FDR) ɑ < 0.5 would be completely meaningless. With 40 tumour and 5 normal samples you should get extremely robust results, and I’d suggest that you should use something like ɑ < 0.0001. Given your PCA and your lack of DE, the only valid conclusion is that you really have no evidence whatsoever for differential expression. The samples are indistinguishable by the state of the art. This in itself is a striking finding.

You can usually get away with FDR < 0.1, but that's as high as you can go. This all presumes you're doing follow-up experiments of some sort, of course. I guess you could increase the FDR more, but you're then really increasing the odds that your follow-up experiments will fail. Obviously increasing the FDR will decrease the confidence in the results, it's called "false discovery rate" for a reason.
Yes, fold-change cut-offs are mostly useful when you need to prioritize hits further. Small changes are almost impossible to validate with less-sensitive methods (qPCR) and less likely to be biologically relevant.
You're just going to waste your time if you ignore FDR.

Remember that all thresholds are intended to increase the odds of success in follow-up experiments. If you have infinite time and resources then you can set these to whatever you want, but you have neither so you better prioritize reasonably.
It's likely that you have a large batch effect in your data, effect sizes in cancer tend to be much larger than this.

Thanks Devon. Yes there will be follow up experiments after selection. So, what would recommend me now? Can I go with FC and FDR < 0.1? Yes, there is batch effect what can I do now? Can I apply subsampling procedure to get more DEGs?

Yes, FDR<0.1 and ideally some reasonable fold-change (say a 30% change). Have a look at packages like sva for dealing with your batch effect.

Is it a good idea to apply subsampling procedure for DEA?

@beginner No, that's not useful, presuming you mean subsampling from the counts.

selecting random samples from cancer cohort and repeating differential analysis n times will give different set of DEGs and then select which DEGs are most repeated in all the analyses. This is the one I'm asking.

No, that's not helpful, you're just approximating the results you would have gotten otherwise.

Sure thank you. And one last question. After selecting candidates..for experiments I have selected two cell-lines.To look which gene is more expressed in which cell line, doing a simple t-test between the cell-lines for those genes is fine or should I go with Deseq2 to check how those genes were expressed?

Use DESeq2, it's more powerful than a T-test.

Yes, it is more powerful but Deseq2 needs replicates for performing differential analysis. It will give a warning/error if tried to analyze data without replicates.

I assumed you were talking about cell lines with replicates. No method will reliably work without replicates.

No, I just have two breast cell-lines. Interested in doing experiments using one cell-line. So I thought of looking at candidates expression in those cell-lines and then select based on that.

@beginner It’s the other way round: DESeq needs fewer replicates than a standard t-test, because it pools variance estimates. This is, in fact, the whole reason why dedicated DE analysis packages exist. With 45 samples a simple t-test should actually work just as well as DESeq though.

@KonradRudolph you have mistaken. please check my previous comments.

@beginner I’m puzzled. First off, your question mentions 45 replicates, not two (or one). Secondly, how do you suggest performing a t-test on two replicates?! (One replicate doesn’t even make sense … what would you compare against?)

No. After doing differential analysis between 40 tumor and 5 normal I get some differential expressed genes. From that I will select strong candidates and do some experiments. For the experiments I have two cell-lines. To select the good cell-line where I can see high expression of those candidates I have to do differential analysis between those two cell-lines. Note: Please have a look at my comments and Devon's comments.

I think since your PCA shows that your normal samples cluster right alongside the tumor ones, you have to conclude that there is much more variation between tumors than between tumors and controls.  DESeq's differentially expressed gene list was generated properly.
Your results are what they are, sorry they are so useless.

In general, messing with FDR or log fold change cutoffs beyond the ones you used gets pretty sketchy. I usually only adjust these values in the more conservative direction, although you may be able to get away with FDR < 0.1 if you have a large effect size (which you don't).
You might consider investigating why you don't seem to get any differentially expressed genes before adjusting your cutoffs: PCA on the samples based on their top 1,000 most variably expressed genes would be a good start. Do the tumor and normal transcriptomes cluster apart from each other?
With so few normal samples, it is possible that you are not fully characterizing the variability of that population and it is affecting your ability to identify meaningful differences between normal and tumor.

Hi, I made the PCA and it is given in the question. please have a look.

If you have a large effect size you wouldn’t need to increase the FDR to 0.1: a lower FDR would still capture the effect.


text:
Merge / Reconciliate several de novo transcriptome assemblies with different kmers
I am building a De Novo transcriptome reference assembly for an eukaryotic organism for which I have a genome.
I've created several assemblies with rnaSpades using different kmer sizes (19 to 69 with step 10). Now I would like to merge them into one final transcriptome.
How could I do that?
Is using a genome assembly reconciliation tool such as metassembler a good idea?

Do you have a genome for the species you have data for?

@conchoecia I do have a genome, although it's only a draft genome. I'm currently using it to compare the quality of my assemblies

Possible? Sure, but there is a danger in doing meta-assembling. BUSCO is an ever-increasing score - if you merge multiple assemblies it will always improve because it's not really telling you much about false positives. I am not sure what do you want to use the transcriptome assembly for, but those won't ever be perfect and loads of effort might lead to only a marginaly better transciptome.

You may combine all your transcriptomes into a single file and then apply a clustering method to group very similar transcripts into a single one. For that purpose, you may try CD-HIT-EST or MMseqs2. For each identity threshold you are going to test, you may assess the final quality with BUSCO or by blasting against reference sequences.

A simple option would be to simply supply SPAdes with some contigs that you like with the --trusted-contigs option.
There exist tools specifically addressing the problem of merging transcriptome assemblies. I am unsure whether this is notably different from genome assemblies sufficiently that they are better than the genome assembly mergers. Here is a partial list:

transfuse
this paper
this paper
this paper
DRAP

This paper seems to have some more information on comparing tools, though it is not focused specifically on merging.
See also the SeqAnswers on this topic, and this discussion too.


text:
Bioinformatics approach to dentifing potential PCR primer sequences for transcribed gene
I have an annotated transcriptome and would like to develop PCR primers for particular transcribed genes. My species is a non-model plant. Can I use BLAST or another tool to identify potential PCR primer sequences?  Or more generally, is there a workable bioinformatics approach to identifying potential PCR primers for transcripts arising from a non-model organism? Any literature refs would be great too.

Aren't classic primer design tools such as https://www.ncbi.nlm.nih.gov/tools/primer-blast/ enough? Is there something special here I am missing?

I guess there are limitations when designing primers within smaller transcripts.

It depends if you want to amplify your target from genomic DNA or cDNA.
Depending on intron size, primers that anneal to different exons may not be able to amplify across the intron(s), and if you're unlucky enough to design a primer that spans a splice junction, you'll definitely be stuffed if you're working with genomic DNA. I would try blastn or translate to an amino acid sequence and use tblastn to see if you can identify the intron/exon structure in a related species. If you don't have a genome sequence available for anything close, I think your only option is to target as small an amplicon as possible and hope for the best.
If you are working with cDNA, any standard primer design tool should do

So taking the blastn approach, one would have to blast against reference genomes (of vascular plants in this case) to try to get the intron structure? Or would you blast on something else?

I'd be tempted to just blast against the full nr database and see what comes up. Unless you know specifically what the closest reference is


text:
TRAL does not find "Phobos result file"
I want to use TRAL to annotate tandem repeats in the reference genome of Caenorhabditis elegans. For this, I need to install some external software, such as Phobos. I've downloaded Phobos and I am using it by typing the path to it's executable like this:
input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line genome.fa output.txt
This works quite well.
In order for TRAL to find Phobos, I changed TRAL's configuration file as explained here.
I've replaced PHOBOS = phobos with PHOBOS = /Users/user/Documents/.../project/input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line
Then, I followed the tutorial Run and parse de novo repeat detection software and typed the following commands in Python 3 in the Terminal:
import os
from tral.sequence import sequence
from tral.paths import PACKAGE_DIRECTORY

proteome_HIV = os.path.join(PACKAGE_DIRECTORY, "examples", "data", "HIV-1_388796.faa")
sequences_HIV = sequence.Sequence.create(file = proteome_HIV, input_format = 'fasta')

So far so good. But the following command creates an error:
tandem_repeats = sequences_HIV[0].detect(denovo = True, detection = {"detectors": ["Phobos"]})

Did not find Phobos result file in /var/folders/z2/rt6qxntd3tg5bjmqlqnm0mdm0000gn/T/tmpvc20y3i4/001/PHOBOS/phobos.o

Note that if I use PHOBOS instead of Phobos in the command, I get the following error:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/sequence.py", line 162, in detect
    **kwargs['detection'])[0]
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 962, in run_detector
    Detectors(detectors, sequence_type)
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 872, in Detectors
    FINDER_LIST.keys()))
Exception: Unknown TR detector supplied (Supplied: ['PHOBOS']. Known TR detectors: dict_keys(['HHrepID', 'Phobos', 'TRED', 'T-REKS', 'TRF', 'TRUST', 'XSTREAM']))

I am working on Mac OS X. I installed TRAL with pip install tral.

I can reproduce this, so I think it might be a TRAL bug. See https://github.com/acg-team/tral/issues/12

Are you running on DNA or protein? In general, TRAL support for DNA is weaker than we'd like it. If you're searching exons, I might try performing the translation first and running in AA mode.

Running on DNA. I am interested in DNA repeats. Mini-and microsatellites are unlikely to code for anything. If I run it on the translation, then the repeats are not bound to appear at the DNA level because of the genetic code degeneracy.

OK. Just make sure you review the results carefully when using TRAL on DNA, as the parameters were trained for AA. We have a student working on improving DNA repeats now, so this should get easier.

TRAL 0.3.5 has bugs running and parsing Phobos. You can follow progress towards a fix on this issue: https://github.com/acg-team/tral/issues/16
Until the issue is fixed, you should remove Phobos from the repeat detector list in your tral config file (~/.tral/config.ini) and use a different repeat detector.
[sequence]
    [[repeat_detection]]
        DNA = TRED, T-REKS, TRF, XSTREAM


text:
GMAJ error: Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
I have installed GMAJ using:
wget http://globin.bx.psu.edu/dist/gmaj/gmaj.zip  -O output/software/gmaj.zip
unzip -u output/software/gmaj.zip -d output/software/gmaj/
rm output/software/gmaj.zip

I ran lastz on two chromosomes from two different assemblies of C. elegans:
lastz c_elegans.PRJNA13758.WS263.genomicI.fa I.fa --format=maf > MAF.maf
and got a MAF file.
Then, I tried to visualise the file using:
java -jar output/software/gmaj/gmaj.jar MAF.maf
but I got the following error:
Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
    at edu.psu.bx.gmaj.MajGui.setDefaults(MajGui.java:272)
    at edu.psu.bx.gmaj.Maj.<init>(Maj.java:53)
    at edu.psu.bx.gmaj.MajMain.main(MajMain.java:86)

Can someone help me resolve this error or suggest an alternative way to visualise my alignment?
I work on Mac OS X.
Cheers.

It's a programming error that should look like this:
StringBuffer sb = (...)
(...)
String s = (String)sb;

while it should be
StringBuffer sb = (...)
(...)
String s = sb.toString();

( see also https://www.biostars.org/p/177992/ ) . You cannot do anything about it but asking the authors of the software.

Sent an email to Cathy Riemer asking for help.


text:
Psipred installation
I am trying to install psipred in a unix server.
Psipred is installed, but fails because it can not find the correct database files.  I have downloaded uniprod90.fasta. Psipred also loads the old legacy ncbi BLAST-2.2.26 and BLAST+.
runpsipred and runpsipredplus have been patched to call the correct version of BLAST.
I have run a few commands to index the uniprod90 database:
pfilt uniref90.fasta > uniref90filt
formatdb -t uniref90filt -i uniref90filt
makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90.fasta

I then ran psipred and got this error:
$ ./psipred/example/example.fasta
Running PSI-BLAST with sequence  ./psipred/example/example.fasta ...
[blastpgp] WARNING: Unable to open uniref90.pin
[blastpgp] WARNING: Unable to open uniref90.pin
FATAL: Error whilst running blastpgp - script terminated!

makeblastdb created 26 pin files, but not uniref90.pin.

Try with:
makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90

The -out is the name of the DB. So if you use uniref90.fasta, it creates uniref90.fasta.pinwhile the script seems to be expectinguniref90.pin`.


text:
Odds ratio calculations in GWAS paper
Posting here because the reddit thread I made didn't gain any traction.
I'm having a little bit of trouble understanding how the odds ratio was calculated in this table, specifically the OR marked in red:

To my understanding, an odds ratio = $$[(\text{# of cases with exposure}) / (\text{# of controls with exposure})] / [(\text{# of cases without exposure}) / (\text{# of controls without exposure})]$$
In the red-highlighted case, my odds ratio is
$$(0.661 / 0.6580) / [(1-0.661) / (1-0.6580)] = 1.0134$$
which is very off from $0.98$. Also, I know I'm using frequencies rather than absolute counts, but it shouldn't matter because if I were to multiply each frequency by the corresponding population values, all the population values would cancel out.
So, what have I done wrong?
Link to the paper

Those odds ratio aren't calculated from the raw allele frequencies. They're calculated from coefficients estimated from complex linear models with a bunch of covariates that are run on subsets of the data then meta-analysed together.
It is a bit surprising that the odds ratio is less than one when the frequency is higher in cases, but the p-value indicates that the OR is not significantly different from 1 so it may be down to statistical noise. It's possible that within populations one allele is more common in cases than in controls but that cases are primarily sampled from populations where that allele is uncommon, causing it to be more common in controls overall.  See here. Alternatively, 1/0.98 is pretty close to the value you calculated, so it's possible that the reference allele got switched around at some point.


text:
size of the pathways for analysis and filtration
I have recently started working on a substance's effect on a cell line in different dosages. for this, there is a tool called bmdexpress2 that I am using. Its input is the normalized counts from RNASeq for each dosage as a big matrix. When it comes to the pathway analysis step, unlike hallmarks of GSEA, this tool uses some databases for defining pathways which involve pathways of even 2-3 genes.
So what I want to discuss here is the strength of small pathways and filtering thresholds. How can we bioinformatically decide on the importance of a pathway made of 3 genes? Should we just filter them out as there are much more comprehensive pathways? Or do their GO Levels matter? Also there are many cases that out of those 3 genes, 1 of them is differentially expressed; is having 1 gene diff. expressed out of X genes in a pathway enough to say that the pathway is enriched?
Here is an example GO term for a small pathway.
And here is an example of the output I get from the tool. all genes platform column is the number of genes involved in that pathway, input genes column is how many of the genes in that pathway are significantly differentiated. I believe the rest is self explanatory.

Looking at the table, can we really say, ESCRT III Complex, for example is a good pathway to focus while having 1 gene out of 13 in a pathway but that change is significant?
Thanks

Why is the GO term relevant? Is that one of the terms of one of the proteins in your pathway? Is it a term that was reported enriched for that pathway? What sort of pathway analysis is this? Does the tool report over/under represented BP GO terms? In what? Does it return a list of pathways with over/under represented diff expressed genes? Presumably, the statistical test will only show something if it is significant, which will be hard with very small pathways, but please [edit] your question and give us some more details about the analysis.

statistically. the tool only returns fisher exact two tail test value and I already filtered my list of pathways according to that. the tool also returns pathways as go terms and the example go term is one of the significantly enriched

There's nothing special about the term. In order to help you, we would need to see the relevant pathway and the numbers involved. As I said before, please [edit] your question and give us more details.

Are you talking about pathway or about genes related to a GO term? It is not the same as pathway, although if you give us your definition of pathway it might be that it fits... If you don't take into account the structure of GO terms you'll get many false positives. Also, could you [edit] your question and include the data in text format? It will be easier to search more information about the genes and GO terms.. Also did you correct by multiple comparisons (as an answer suggest)?

The problem here is that any one GO term with a small number of members is unlikely to include any differentially expressed genes by chance, but if you test a large number of GO terms, some of them definitley will. The solution is to correct the p-values for multiple comparisons. If you do this, then none of those GO terms with a single DE member will be significant. The larger GO terms with multiple DE members in the chart you show probably won't survive correction either. Hopefully you have some that are more significant.


text:
What is the current state-of-the-art in assembling hybrid transcriptomes?
We are considering attempting de novo assembly of a species transcriptomes (i.e. without a reference genome) using the combined NGS outputs of Iso-seq and Illumina.
One example I saw (Li et al 2017), used the standard PacBio tools to assemble a transcriptome, followed by correcting the assembled sequence using a tool called proovread, followed by cd-hit-est and then cogent.
Would this be considered a reasonable pipeline or is there an alterative recommendation?

My recommendation would be to run Trinity using the --long_reads option, which allows you to provide an error-corrected fastq file for anchoring reads to a transcript isoform. The reads are then clustered per-isoform, and assembled, in a similar fashion to what's done with genome-guided Trinity:
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity
Trinity --seqType fq --max_memory 50G --long_reads corrected_pacbio_reads.fa --left reads_1.fq  --right reads_2.fq --CPU 6

Disclaimer: I have contributed code to Trinity at some time in the distant past


text:
Convert Cotton Probe ID to Gene Symbol
I am new to bioinformatics, my background is in Electrical Engineering. I am trying to convert Affymetrix Cotton Probe IDs to gene symbols. I have a gene expression dataset and I need the expressions for only certain genes. So in the dataset, I have a probe id such as Ghi.10301.1.S1_s_at, how can I get the gene symbol associated with it? Thanks
EDIT: This is the specific dataset I am looking at
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50770

So which file are you trying to parse? Is it GSE50770_series_matrix.txt? The link you gave has various resources associated with it .

One problem you may face is with the state of the cotton genome (Gossypium hirsutum) itself, there are no gene symbols for this genome, you're stuck with accession IDs for each probe, which the provided Affymetrix annotation file has some information on via homology to other characterised proteins and coding sequences, that's about it!  The official genome annotation is not much use either!

Your dataset was produced on a microarray. The probe IDs in the FASTA file are similar to read IDs produced from other NGS experiments like RNA-Seq. According to the workflow provided by Ensembl the first step is to align your reads to a reference. You'll only keep the reads that match the reference completely. Once the probes are aligned you'll be able to quantify the number of reads mapping to the gene that the probe was designed to attract.
I did find a file providing the database entries associated with each probe ID (one of which is the ID you provided) for the Affymetrix Cotton Genome Array.

Hi @kohlkopf thanks for your response. I don't need to quantify the reads here, since I already have expression values for different genes in the dataset. What I am trying to do is convert the Probe IDs to gene symbol so I can extract the data for a specific subset of genes I am interested in. I will add a link of the data set I am looking at on the original post.

This is demonstrably not the answer to the OPs question, the OP has not got a FASTA file of probe sequences, equally the hits < 100 step and "quantify the number of reads mapping to the gene" - would only be appropriate for a RNASeq experiment, in this case expression is via probe intensity measurements not read count, counting the number of probes would simply be a bad proxy for gene length.  This is simply an ID mapping issue.

Mapping to Ensembl won't help either as it only has the Gossypium raimondii genome were as the Affy Cotton array is based on Gossypium hirsutum.

The first stop for decoding probe IDs for a microarray platform is usually the GEO Platform Description: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL8672
This description contains a table of annotation, and is provided (in this instance) by the array manufacturer. You can view, or download, the full table: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL8672&id=1538&db=GeoDb_blob35
GEO Platform annotations do date, and so for up-to-date annotations, download from NetAffx: http://www.affymetrix.com/support/technical/byproduct.affx?product=cotton
The most recent annotations (dated 2016-04-13) are available here: http://www.affymetrix.com/Auth/analysis/downloads/na36/ivt/Cotton.na36.annot.csv.zip.
For the example in your question, looking at this probe in the CSV:
"Ghi.10301.1.S1_s_at","Cotton Genome Array","Gossypium hirsutum","Mar 30, 2016","Consensus sequence","GenBank","Ghi.10301.1","gb:DT050231 /DB_XREF=gi:72275634 /DB_XREF=COT_DW_H07 /TID=Ghi.10301.1 /CNT=21 /FEA=EST /TIER=ConsEnd /STK=1 /UG=Ghi.10301 /UG_TITLE=Transcribed locus, moderately similar to NP_568148.1 ubiquitin conjugating enzyme ubiquitin-like activating enzyme (Arabidopsis thaliana) /REP_ORG=G. hirsutum","DT050231","Ghi.10301","Ghi.10301","---","---","---","---","---","est","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","This probe set was annotated using the Accession mapped clusters based pipeline to a UniGene identifier using 21 transcripts. // false // Accession mapped clusters // E","AI730871,DN803875,DR452973,DR457278,DT050231,DT545642,DT546719,DT546982,DT546983,DT547780,DT547781,DT555433,DT559909,DT564450,DT566924,DT569293,DT569974,DT571307,DT571318,DT571472,DT573548","AI730871 // /DB_XREF=BNLGHi8130 // gb // --- // --- /// DN803875 // /DB_XREF=G.hir-8-10 DAA bolls drought stressed751 // gb // --- // --- /// DR452973 // /DB_XREF=CM008C10 // gb // --- // --- /// DR457278 // /DB_XREF=CM053F01 // gb // --- // --- /// DT050231 // /DB_XREF=COT_DW_H07 // gb // --- // --- /// DT545642 // /DB_XREF=EST1056282 // gb // --- // --- /// DT546719 // /DB_XREF=EST1057359 // gb // --- // --- /// DT546982 // /DB_XREF=EST1057622 // gb // --- // --- /// DT546983 // /DB_XREF=EST1057623 // gb // --- // --- /// DT547780 // /DB_XREF=EST1058420 // gb // --- // --- /// DT547781 // /DB_XREF=EST1058421 // gb // --- // --- /// DT555433 // /DB_XREF=EST1066073 // gb // --- // --- /// DT559909 // /DB_XREF=EST1070549 // gb // --- // --- /// DT564450 // /DB_XREF=EST1075090 // gb // --- // --- /// DT566924 // /DB_XREF=EST1077564 // gb // --- // --- /// DT569293 // /DB_XREF=EST1079933 // gb // --- // --- /// DT569974 // /DB_XREF=EST1080614 // gb // --- // --- /// DT571307 // /DB_XREF=EST1081947 // gb // --- // --- /// DT571318 // /DB_XREF=EST1081958 // gb // --- // --- /// DT571472 // /DB_XREF=EST1082112 // gb // --- // --- /// DT573548 // /DB_XREF=EST1084188 // gb // --- // ---","---"


text:
Effect of mutation in DNA sequence on transcription factor binding sites
How much does a single mutation/alteration of a nucleotide affect the presence of a transcription factor binding site (TFBS)?
I am from computer science background(Obviously).
I want to make a general assumption about the number of mutated bases in a TFBS to minimize the binding affinity?

I haven't worked on TFBSs but I really, really doubt you will be able to do this. In mutations, like in real estate, it's location, location, location. It isn't a question of how many mutations but what mutations and where. TFBSs are usually described using position weight matrices, precisely because that's how you can show the importance of different changes in different positions.
Here are a couple of random Google image search results for "transcription factor binding sites logo":

(This was taken from here and you might want to read through it.)

(source: Front. Genet., 23 February 2016 | https://doi.org/10.3389/fgene.2016.00024)
The main take home messages from these images are:

Different positions have different importance and potential to disrupt binding. In the first image, I would expect that changing the T at position -3 to anything else would be pretty bad. However, changing an A to a T at -8 would likely be neutral.
Different TFs have very different binding sites and very different binding site lengths.  Look at panel B of the second image. That one appears to be a very conservative site and likely any change would be deleterious. Conversely, changing any of the nucleotides of the 17bp spacer shown in panel A would likely not be important as long as you don't change the length.

Please note that I have not done more than glance at the sources I took these from, so my inferences may well be wrong, I only use them to illustrate the general point: you cannot really put a number to this. You will always need to look at each specific change and each specific position. In some cases, a single change can disrupt binding while in others multiple changes might have no effect at all.

You might want to look into TFBS databases like JASPAR or TRANSFAC.

Every motif has a score. The score is calculated from the number of times each base was observed in that motif in the genome (I'm not sure how exactly these are calculated). Higher numbers indicate that this motif was found more frequently, which we take as a proxy for it providing better binding.
The Ensembl VEP can give you the score change for a mutated motif in the output, so this may be the easiest way to analyse this.

The motifbreakR package for R I believe does something like what you want to do;
https://www.rdocumentation.org/packages/motifbreakR/versions/1.2.2/topics/motifbreakR
If you're looking for a general rule as to what kind of mutations/where in the motif mutations will significantly affect motif binding, I suggest generating a lot of sequences with simulated SNP's (single nucleotide polymorphisms) and seeing how these score against the particular motif. Then, you will have a large(ish) database of sequences with one base pair difference from the motif and their corresponding score. This will allow you to judge what kinds of mutations affect the score the most.
Furthermore, you could somehow make a distribution of scores versus mutations and generate an empirical p-value to determine statistical significance (for instance, such and such mutation caused a score that was in the lowest 5% of scores for all mutations). You would know better than me how to do this, since your background is in CS and my PhD subject is genetics.
To make an assumption about the number of bases that need to be mutated, do this for one, two, three, etc. simulated SNP's.


text:
Merging ssu-align alignments
ssu-align is an HMM-based MSA tool for ribosome's small subunit (SSU) sequences. It's bundled with three SSU models:

archaeal,
bacterial
eukaryotic

Approach Given a set of archaeal and bacterial sequences the tool will first identify the best model for each sequence and produce two separate alignments for archaea and bacteria, and the tool provides no merging functionality.
Rationale The underlying issue with SSU is the secondary structure, which results in: 1) indels; 2) pairwise point mutations (to maintain the secondary structure)
Problem While you can force it to use a single model for all sequences, aligning archaeal sequences to the bacterial model will produce poor results.
Solution? Nevertheless, the authors of GreenGenes claim that they've used ssu-align to produce two MSAs (one for archaea and bacteria) and build a joint phylogenetic tree using FastTree, but they say nothing about the way they've managed to merge the alignments. Am I missing some functionality from SSU-align?
Note Merging two giant MSAs using conventional tools, such as MAFFT, won't get me far.
P.S.
I've already tried contacting the authors of the GreenGenes paper, but the corresponding author hasn't responded.

I think this is more to do with tree building than alignment construction.
I believe the underlying issue is the 1) distribution of indels, this results in a blockwise deletion in the resulting phylogeny and 2) nucleotide saturation.
Indels Thus, for indels
Archaeal
XXXX----YYYYY
bacterial
----ZZZZYYYYY

The resulting phylogeny will be based on YYYYY and XXXX is ignored.
The alternative explanation is that the genetic distances involved have resulted in saturation of the nucleotide sequences. SSU has to be a nucleotide phylogeny and the level of reversion mutations over a vast amount of evolutionary time can become indistinguishable from random.
With regards the first point, you can replace the indels "-" with "N". This is a cheat to trick a tree building program into not discarding the indel. Model based tree building programs have no ability to create a genetic distance between a base pair and an indel, so this data is essentially discarded from the tree. Worse it will remove the site from phylogenetic inference - which can be a massive loss of data. Two notes one this,

Some max. likelihood, Bayesian tree programs can avoid blockwise deletion, I don't know programmatically how they do this;
In PAML the distribution of indels was modelled, it doesn't appear to have been widely adopted however, but might be of value.

Alternative The alternative method is really for protein genes and called "RY" coding, it might work on SSU because of the compensatory mutations involved in maintaining a secondary structure of the SSU.
Word of warning, the saturation may be so great, very little nucleotide point mutation is preserved. Essentially we are talking about preservation of phylogenetic signal across the origins of life ... and that is a very, very long time ago.


text:
How to run MaxQuant in command line mode?
MaxQuant is a software package for mass spectroscopy and proteomics. There is a windows version and a linux version. To run on linux you have to use a program that is called mono. I think, it is developed by Microsoft, which I find quite nice from them. I installed mono on Linux and started MaxQuant.
mono MaxQuant.exe

The documentation claims it is for running MaxQuant in commandline mode, though a GUI is started anyway. How can I start the program in commandline mode to perform batch analysis? I think, there is not enough information given on the web-page.

Give
mono MaxQuantCmd.exe [OPTIONS]

a try.
University of Oslo description of how to use the tool in their infrastructure.

Kohl, the link is broken. Could you update it please?

You can now install MaxQuant simply with Conda:
conda install -c bioconda maxquant

and to run MaxQuant type:
maxquant mqpar.xml

Reference


text:
Using sickle for quality trimming
I am using sickle package from bioconda to trim my reads from the Illumina Miseq.
When I tried following code, it failed to run
sickle pe -t illumina -f file1   -r file2  -o output1  -p output2  -s combined_output -q 27 -l 25

When I replaced illumina with sanger, it works. But my reads are from the illumina. Could anyone clarify if using this option is acceptable ? Also, what is the optimum value of -q and -l

Sickle examines quality and length quality thresholds via a sliding window. This assesses whether to trim the 3'-end of reads, for example removing a sudden drop in quality. It can also assess whether the quality of the 5' end is sufficiently high for trimming.

It'd be rare to run into Illumina encoded quality scores any more, I doubt there are machines running that still produce those. Everything produced at the moment uses Sanger encoded (phred + 33) scores.
There is no optimum value for -q and -l, though reasonable values would typically be <=5 for -q and ~5 for -l.


text:
Simplest way to work out structural variant type?
In VCF 4.2, a structural variant (SV) can be described with the BND keyword in SVTYPE. For example, the following example is an insertion (from https://samtools.github.io/hts-specs/VCFv4.2.pdf):
#CHROM    POS ID REF ALT                      QUAL FILTER INFO
2 321682  bnd V  T   ]13 : 123456]AGTNNNNNCAT 6    PASS   SVTYPE=BND;MATEID=bnd U
13 123456 bnd U  C   CAGTNNNNNCA[2 : 321682[  6    PASS   SVTYPE=BND;MATEID=bnd V

The keyword can be used to describe any SV event.
I would like to work out how to find out if an event is an insertion, deletion, duplication, inversion, etc. The algorithm/program will tell me the example above is an insertion.
I'm not able to find a tool that can do the conversion for me.
PS: The program I'm using is GRIDSS. This program would only generate BND in the SVTYPE field. It's not strictly a bug (BND can be used for any SV arrangement), but it's annoying in analysis.

how about 'just' using the INFO/SVTYPEattribute ?

@Pierre SVTYPE is BND.

@Pierre Some SV programs such as GRIDSS only output BND to the SVTYPE. Otherwise, I would not be here asking.

from the VCF spec: INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant"> Value should be one of DEL, INS, DUP, INV, CNV, BND.. So BND (break end) is NOT an INDEL, and NOT an insertion.

@Pierre My above example took from the specification.

"Some SV programs such as GRIDSS only output BND to the SVTYPE" oh, I see

May I suggest you simply not use GRIDSS and instead use a tool that produces valid output?

@terdon What do you mean? GRIDSS is a published SV program. BND is a valid VCF output?

If it uses BND to indicate an insertion, then it is not following the VCF specs since the SVTYPE for insertions should be INS and BND should only be used for break ends.

GRIDSS's author had posted his own solution on Github. The code is:
simpleEventType <- function(gr) {
    return(ifelse(seqnames(gr) != seqnames(partner(gr)), "ITX", # inter-chromosomosal
           ifelse(gr$insLen >= abs(gr$svLen) * 0.7, "INS",
           ifelse(strand(gr) == strand(partner(gr)), "INV",
           ifelse(xor(start(gr) < start(partner(gr)), strand(gr) == "-"), "DEL", "DUP")))))
}

Can you please elaborate on this answer? What's gr? Looks like R, which library should I use for this?


text:
Get Gene Expression Matrix from GEOquery
I am trying to get gene expression matrix for a list of genes
I have for my list: Gene ID, Gene Symbol
How can I get for each gene corresponding expression as array??
I suppose that I will have
GENE_ID GENE_NAME SAMPLE_1 SAMPLE_2 .....
I am new in that, I don't know more details
I am trying to study relation between diseases and genes(human) using gene expression as features.
Thanks in advance.

Are you working with a particular dataset?

I am using GEOquery library

The answer to this really depends on the type of data you're retrieving from GEO. Microarray data sets should have a normalised matrix of expression values uploaded as part of the entry. getGEO() from the GEOquery package will return a list of ExpressionSets - in many cases this list will have a length of 1. An ExpressionSet has an exprs() accessor for retrieving the expression matrix. An example (GSE76510, an Illumina HT12 v4 BeadArray dataset):
> library(GEOquery)
> gse76510 <- getGEO('GSE76510')
> names(gse76510)
[1] "GSE76510_series_matrix.txt.gz"
> class(gse76510$GSE76510_series_matrix.txt.gz)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
> exprs(gse76510$GSE76510_series_matrix.txt.gz) %>% head
             GSM2027280 GSM2027281 GSM2027282 GSM2027283 GSM2027284 GSM2027285 GSM2027286
ILMN_1343291  14.936287  14.983927  15.008199  14.830549  15.008199  14.789385  14.891532
ILMN_1343295  13.404168  13.329903  13.289727  13.222221  13.386299  13.240755  13.125056
ILMN_1651199   6.779805   7.189868   7.046952   6.962957   6.697924   6.850950   6.792647
ILMN_1651209   6.983803   6.769702   7.016746   6.954751   6.890472   6.982947   6.915671
ILMN_1651210   6.837042   6.999098   6.895225   6.944250   6.955914   6.917959   6.837368
ILMN_1651221   7.034363   7.025626   6.606457   6.977622   7.051120   6.971884   7.066445
             GSM2027287 GSM2027288 GSM2027289 GSM2027290 GSM2027291
ILMN_1343291  14.891532  14.983927  14.848452  14.913474  15.008199
ILMN_1343295  13.240755  13.010757  13.525999  13.009661  12.998406
ILMN_1651199   6.786784   6.820451   6.851160   7.013586   6.864316
ILMN_1651209   6.983701   6.995199   7.270760   6.901377   6.950876
ILMN_1651210   6.804955   7.033074   6.894117   6.770691   7.000786
ILMN_1651221   7.003035   6.995199   7.050040   6.923861   6.929933

You will be able to filter out a list of required genes using the probe annotation (available in the appropriate annotation.db - in this case, for example, illuminaHumanv4.db).
If the dataset you're interested in is RNA-Seq, you're probably out of luck. There's currently no requirement on data submitters to provide expression estimates with an RNA-Seq dataset, so the ExpressionSet that GEOquery returns will most likely have an empty exprs slot.
Again, an example (GSE117133):
> gse117133 <- getGEO('GSE117133')
> exprs(gse117133$GSE117133_series_matrix.txt.gz)
 GSM3271892 GSM3271893 GSM3271894 GSM3271895 GSM3271896 GSM3271897

Sir, I am confused which dataset should I use that contains the list of my interest and  help me in my model to predict relation between genes and diseases. Do you suggest to use GSE76510 (MY list is human proteins that included in DrugBank DB)?? Thanks in advance.

@EmanIsmail I’m not suggesting any specific dataset here - they are merely to illustrate how the functions work in R. You need to find a dataset in GEO that assays gene expression in your system of interest.


text:
rJava issue when trying to use ATACseq pipeline
I am trying to use the following ATACseq pipeline.
Unfortunately when I am importing it with library(esATAC), it is giving me an error:
Error: package or namespace load failed for ‘rJava’:
 .onLoad failed in loadNamespace() for 'rJava', details:
  call: dyn.load(file, DLLpath = DLLpath, ...)
  error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so':
  dlopen(/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so, 6): Library not loaded: /Library/Java/JavaVirtualMachines/jdk-9.jdk/Contents/Home/lib/server/libjvm.dylib
  Referenced from: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so
  Reason: image not found

There are threads that discuss similar problems 1 and 2:
But I checked the Java version and it is the right one: 64bit:
d-172-25-60-75:~ nikitavlasenko$ java -version
java version "1.8.0_172"
Java(TM) SE Runtime Environment (build 1.8.0_172-b11)
Java HotSpot(TM) 64-Bit Server VM (build 25.172-b11, mixed mode)

I tried to do the following in Rstudio:
Sys.setenv(JAVA_HOME = "/usr/bin/java")
library(rJava)
library(esATAC)

But library(rJava) fails. I am using macOS High Sierra 10.13.4. I am not sure what can I try further. Any suggestions would be greatly appreciated.

Fundamentally, JAVA_HOME needs to be set to the Java home, not to the Java binary. This is a folder path that ends in Home.

I uninstalled Java the way it is written here:

Click on the Finder icon located in your dock
Click on Applications tab on the sidebar
In the Search box enter JavaAppletPlugin.plugin
This will find the JavaAppletPlugin.plugin file
Right click on JavaAppletPlugin.plugin and select Move to Trash

Then installed the recommended latest version for mac on the official Java website and it fixed the issue.


text:
Supplying adapter sequences for ATACseq pipeline
I have several atac_seq.gz files and I need to run the pipeline with them.
I have single end ATACseq, so the function atacPipe requires me to have some kind of adapter sequence for it to run. What is it? Where could I get one?

Update

I found the adapter sequences, supplied them, but it is writing me that

Error in atacPipe(case = list("mm10", fastqInput1 = paste0(base_folder,  :
    adapter1 should not be NULL for single end sequencing data

library(esATAC)
base_folder <- "/Users/nikitavlasenko/Downloads/sams_dev/"
conclusion <- atacPipe(case=list("mm10", fastqInput1 = paste0(base_folder, "Micro_atac12.fastq.gz"),
                   adapter1 =
"CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT",
                   createReport = TRUE))

I supplied it and it is not null, so I am not sure what I am doing wrong.

Presumably one would have to know which platform the data came from since (at least for RNA seq) there are different adapter sequence sets. This write up has an adapter sequence agnostic approach.

I was calling the function incorrectly. Here is how it worked:
conclusion <- atacPipe("mm10", fastqInput1 = paste0(base_folder, "Micro_atac12.fastq.gz"),
                   adapter1 = "CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT",
                   createReport = TRUE)


text:
Genome scaffolding
I have assembled a virus genome using Ray resulting in approximately 5000 contigs. Now I want to build a scaffold of those contigs to get the full genome (I am aware of the fact that there might be some gaps in the genome).
I  found one tool called Medusa to do this. Unfortunately, this tool cannot take a file above 50 MB and my contig file is about 156 MB.
Are there any other tools that can perform reference based scaffolding?
It would be great if anyone knows about genome scaffolding packages in Bioconda.

Any genomics tool that doesn't work over 50MB data is meaningless. Please delete Medusa.

@SmallChess if it's designed for viral genomes, why not.

The number of contigs and total assembly size you have suggest that there was probably more in the sequencing run than a single virus strain. Does the total assembly length correspond to expected genome size of the reference strain? Have you tried just to blast some your contigs? Could you have a contamination or something?
Maybe you can just sort out the contigs by the reference using something like artemis or quast compares your assembly to reference and showing nice stats (just to be sure that your assembly makes sense so far). Or you could use directly MUMmer (which the alignment tool that is running on the backend of of Medusa and Quast).
This is not directly answering your question "How to scaffold", however I think you should first figure out if your contigs make sense before you scaffold them. Also sorting your contigs by artemis using reference is practically the same thing as scaffolding since you won't know the sequence in between in either of the cases.

You are right. There is a lot of host genome contamination. I blasted all the contigs against reference genome and only about 10 contigs had significant match with the reference. I was looking for some tool that could take these 10 contigs, re-order them and join them together or add NNNNs if the sequence is missing between the contigs.

And do you know the host genome? If yes, you could map everything to both host and the reference and then kick out everything that maps to host and scaffold just non-host contigs.

I know the host genome because it is cell-cultured virus. And I am actually looking for some tools that could actually do what you have mentioned.

I would personally just run MUMmer of the assembly vs catinated reference of the host and related strain. MUMmer allows secondary alignments as well so you will be able to find also regions that are similar between the both, then parse the output in python/R/whateveryoulike to create a list of scaffolds you want keep/remove and use it to filter out our reference (a script to do that could look like this https://github.com/KamilSJaron/generic_genomics/blob/master/fasta2extract_by_list_of_headers.py). I should probably add this to the answer.

Thanks Kamil. I will try that script. I also found that Viral NGS Pipeline from broad institute can actually do this. There is a python script called order_and_orient which serves this purpose. Here is the link https://viral-ngs.readthedocs.io/en/latest/assembly.html

I found that Viral NGS Pipeline from broad institute can actually do this. There is a python script called order_and_orient which serves this purpose. Here is the link  viral ngs assembly

What you mean is assembly based scaffolding, as opposed to using reads with long distance information such as mate pair/long jumping distance to scaffold.
You could reduce your contig file to just the longest contigs, since they have the most information. Probably 156MB of data is overkill for your virus.
There are other programs out there.
https://github.com/ksahlin/BESST
https://github.com/institut-de-genomique/MaGuS
Your mileage will vary.

I tried that. I took five largest contigs. The output is not good. I will try the tools that you provided. Thank you.

Isn't BESST only for read-based scaffolding?