text:
What process and input data is required for a cellranger reference transcriptome?
I'm analysing single-cell RNA-Seq data using the 10X Genomics cellranger platform. While they provide reference data for Mouse and Humans, other species require a reference "transcriptome" to be generated with the cellranger mkref tool. However, I've found it to be very temperamental giving errors such as:
gffread genome_GFF3_genes.gff -T -o genome_genes.gtf

cellranger mkgtf genome_genes.gtf genome_genes.filtered.gtf \
                   --attribute=gene_biotype:protein_coding \
                   --attribute=gene_biotype:lincRNA \
                   --attribute=gene_biotype:antisense \
                   --attribute=gene_biotype:IG_LV_gene \
                   --attribute=gene_biotype:IG_V_gene \
                   --attribute=gene_biotype:IG_V_pseudogene \
                   --attribute=gene_biotype:IG_D_gene \
                   --attribute=gene_biotype:IG_J_gene \
                   --attribute=gene_biotype:IG_J_pseudogene \
                   --attribute=gene_biotype:IG_C_gene \
                   --attribute=gene_biotype:IG_C_pseudogene \
                   --attribute=gene_biotype:TR_V_gene \
                   --attribute=gene_biotype:TR_V_pseudogene \
                   --attribute=gene_biotype:TR_D_gene \
                   --attribute=gene_biotype:TR_J_gene \
                   --attribute=gene_biotype:TR_J_pseudogene \
                   --attribute=gene_biotype:TR_C_gene

cellranger mkref --genome="reference-name" \
                 --fasta="transcriptome.fasta" \
                 --genes="genome_genes.filtered.gtf"

Fatal INPUT FILE error, no valid exon lines in the GTF file: reference_name/genes/genes.gtf
  Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA file.
subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status 104
subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status -11

What input data is required for the cellranger to construct a reference? What format do the fasta and gtf files need be in? Is it possible to use a "transcriptome" rather than a "genome" as a reference here? For example, I have a transcriptome reference with fasta names >name for each transcript rather than each chromosome.

Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with genome_genes.filtered.gtf, it needs the genomic sequence to go along with that.
As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look at the gene_biotypes in your file and decide for yourself what (if anything) needs to be filtered. For example, I made a cellranger reference last night for dm3 (the old fruit fly genome) with:
cellranger mkgtf annotation.gtf genes.gtf \
    --attribute=gene_biotype:protein_coding \
    --attribute=gene_biotype:ncRNA \
    --attribute=gene_biotype:snRNA \
    --attribute=gene_biotype:snoRNA \
    --attribute=gene_biotype:pre_miRNA

That should be customized for relevance. It's quite possible that your GTF doesn't even have gene_biotype attributes, that's specific to Ensembl I think.

Indeed, I am not working with Humans on this project (I tried to include a "model-organisms" tag to emphasise this). I included this from the cellranger documentation while troubleshooting the issue. To clarify, the reference "genome" is a FASTA with chromosomes named? I have a recent annotation (GTF) for A. thaliana but the reference (FASTA) is by transcripts in this database. An older full-genome reference is available but I'd rather use a recently updated on is possible.

Yes, the reference genome is the fasta file containing chromosomal entries. You can get appropriate GTF files from Ensembl, a transcriptome fasta file is of no use to you.


text:
How to normalise scRNASeq data for differential expression analysis
I wish to perform differential expression analysis for cluster-specific gene expression in single-cell data (with a tool such as MAST or SCDE).
I have data for 3 biological replicates. I performed analysis on 10x data with cellranger and have filtered out cells with low UMIs. I have corrected for batch effects and identified clusters by performing graph-based clustering on the UMAP projections.
As I understand it, the components from batch effect correction algorithms are only intended to be used for the purpose of plotting and clustering. They are not to be used for analysis of individual gene expression or differential expression. This has been stated unequivocally by the original authors of MNN and CCA batch correction algorithms.
I have derived a UMI count matrix for each sample which corrects for amplification bias. Differences between the cells and samples still persist in the raw data. These are obvious from the mean UMI counts to the UMAP of uncorrected data. The number of UMIs differs considerably between samples (batches) and individual cells within a sample (the total UMI per cell). The amount of data from each of the cells varies considerably. The read-depth between samples is also different. Therefore, I need a simple approach to correct for differences in read-depth, total UMI or sensitivity between cells.
Is it possible to normalise UMI counts between cells? Should UMI count data be treated differently to reads per transcript in bulk RNA-Seq? Would it be more appropriate to do so at the level of individual cells or for batches across samples?

I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do:
markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch')
(change "object" and "batch" accordingly)
See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's group for some more info on the LR test.

Thanks, Tim. I was considering to use Seurat. Do I need to normalise the data before doing this or does the ScaleData function take care of that for you?

You will need to log-normalize the UMI counts, you don't need to run ScaleData if you've already found your clusters using the batch-corrected data

source: Stack Exchange

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

sort: votes


text:
Can open-source software be peer-reviewed and published?
My colleague and I have developed a software tool intend to release it open-source. This tool is specifically for tasks in bioinformatics but we think it would be helpful for the wider community. Our institution will permit us to release it provided we get appropriate credit.
Thus we wish to publish it in peer-review. Is peer-review publication of bioinformatics software available? If so what is required to publish it?

There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck publishing in the Nucleic Acids Research yearly special webservers issue, since we make Galaxy wrappers around our tools.
You can also publish on tools in venues like Nature Methods, but it's vastly more difficult to get in there if your paper is purely on a tool and it's not something game-changing like salmon or kallisto.
Regarding what's actually needed to publish a tool, you will typically need the following:

Decent documentation. I ask for this as a reviewer and am asked about it when I submit things.
Easy installation. No one wants to spend two hours trying to compile your tool, so make sure there's a conda package or at least a docker container for it.
Test data. Without this your reviewers and users will likely not be happy since they won't be able to try it out on something small.

For the actual submission you'll tend to need comparisons to other tools. I'm generally not a big fan of this, but you will find that some reviewers will demand it (assuming there are other tools that do something similarish to yours).

Yes, there is even a journal dedicated to this: The Journal of Open Source Software
Here the software itself is peer reviewed and a decent documentation and tests are required. Since the submission and review process is based on GitHub, publishing there is straightforward if your code is on GitHub already.
Note: I am not affiliated with this journal, I just published two software packages there a while ago.

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


text:
What is cellranger doing in comparison to other methods?
I've recently started working with the 10X-Genomics platform with Illumina (MiSeq and HiSeq) for single-cell RNA-Seq. I've been recommended the "cellranger" (version 2.1.0) which I understand handles the barcoding of the platform and performs processing of the data. For example:
1) convert raw Illumina .bcl files into .fastq (demultiplex)
cellranger mkfastq

2) perform alignment, barcode counting, and generate gene-barcode-martices
cellranger count

However, it is unclear which underlying bioinformatics procedure (with which parameters) are being used here. For this pipeline, I wish to know how the trimming, alignment and counts are being performed. This is important to understand whether further normalisation and quality checking is necessary. What steps would I need to do to perform a similar analysis?
If for example, I wish to document in more detail what has been performed here in a reproducible manner or reproduce the analysis pipeline on a different platform (such as Nadia) to compare them.
Ideally, I would like to know the analogous pipeline using open-source tools but at the moment I am first seeking to understand what is going on under the hood when I am running the cellranger.

Have you looked at the cellranger source code (e.g., the alignment command is generated here)? If you want the exact commands run that's the easiest method.

Thanks for the link to the source code. That is helpful but full for the benefit of those not well-versed in Python (or these tools) more explanation of what these parameters are doing would be helpful.

As an update: kallisto | bustools, dropEst, and DropSeqPipe are valid alternatives for this task.

cellranger mkfastq is not necessary anymore.  It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself.  Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.
Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene.  You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam
If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses).  I have put the same data sets through both 10X and Dropseq, and they return about the same thing.

Great answer but it would help to have more detail on how the cellranger tool is performing the alignment, is this with the default parameters for STAR and samtools? I've been advised to "demultiplex" MiSeq output but that may be for older versions of cellranger.

The best approach would be to contact 10x Genomics support: support@10xgenomics.com . They usually respond within a few hours.
Otherwise, you can check the methods from the first 10x Genomics paper:

The Cell Ranger Single-Cell Software Suite was used to perform sample
  demultiplexing, barcode processing and single-cell 3′ gene counting
  (http://software.10xgenomics.com/single-cell/overview/welcome). First,
  sample demultiplexing was performed based on the 8 bp sample index
  read to generate FASTQs for the Read1 and Read2 paired-end reads, as
  well as the 14 bp GemCode barcode. Ten basepair UMI tags were
  extracted from Read2 (14 libraries were made with 5 bp UMI tags, as
  noted in Supplementary Table 1, due to an earlier iteration of the
  methods. For these samples, 5 bp UMI tags were extracted from Read2.).
  Then, Read1, which contains the cDNA insert, was aligned to an
  appropriate reference genome using STAR35. For mouse cells, mm10 was
  used and for human cells, hg19 was used. For samples with mouse and
  human cell mixtures, the union of hg19 and mm10 were used. For ERCC
  samples, ERCC reference
  (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt)
  was used.
Next, GemCode barcodes and UMIs were filtered. All of the known listed
  of barcodes that are 1-Hamming-distance away from an observed barcode
  are considered. Then, the posterior probability that the observed
  barcode was produced by a sequencing error is computed, given the base
  qualities of the observed barcode and the prior probability of
  observing the candidate barcode (taken from the overall barcode count
  distribution). If the posterior probability for any candidate barcode
  is at least 0.975, then the barcode is corrected to the candidate
  barcode with the highest posterior probability. If all candidate
  sequences are equally probable, then the one appearing first by
  lexical order is picked.
UMIs with sequencing quality score >10 were considered valid if they
  were not homopolymers. Qual=10 implies 90% base call accuracy. A UMI
  that is 1-Hamming-distance away from another UMI (with more reads) for
  the same cell barcode and gene is corrected to the UMI with more
  reads. This approach is nearly identical to that in Jaitin et al.4,
  and is similar to that in Klein et al.8 (although Klein et al.8 also
  used UMIs to resolve multimapped reads, which was not implemented
  here).
Last, PCR duplicates were marked if two sets of read pairs shared a
  barcode sequence, a UMI tag, and a gene ID (Ensembl GTFs GRCh37.82,
  ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz
  and GRCm38.84,
  ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz,
  were used). Only confidently mapped (MAPQ=255), non-PCR duplicates
  with valid barcodes and UMIs were used to generate gene-barcode
  matrix.
Cell barcodes were determined based on distribution of UMI counts. All
  top barcodes within the same order of magnitude (>10% of the top nth
  barcode, where n is 1% of the expected recovered cell count) were
  considered cell barcodes. Number of reads that provide meaningful
  information is calculated as the product of four metrics: (1) valid
  barcodes; (2) valid UMI; (3) associated with a cell barcode; and (4)
  confidently mapped to exons.

Of course, they have since updated the library prep and the Cell Ranger software, so the exact protocol is somewhat different.


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

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

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

OUTPUT
TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427
However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?
I'm left with the output however I've no idea how to link this with possibly an online source like NCBI to actually find the source of my ORF (i.e. which protein it comes from). Is there anyway to use Biopython for this as well?

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


text:
Proteins that cannot form biofilm?
I am trying to build a machine learning training set for bacterial protein sequences that form biofilm, and that cannot. I collected the positive sequences from the GO ontology website but for negative sequences I am not sure which sequences to incorporate into my training set.
Can anyone point me to resources for proteins sequences that are known to be not capable of forming biofilms?

Maybe I am wrong but all the proteins you don't have as your positive set is your negative set. How do you select the biofilm capable proteins? By those annotated with a specific GO term?

One option is to treat all proteins that are not explicitly known to form a biofilm as proteins that are incapable of doing so, but this will likely result in many false negatives.
You could take a look at one-class classification and positive-unlabeled (PU) learning. Both are techniques specifically designed for these kind of problems where only one of the two classes is known. (In this case: only biofilm forming proteins are known.)

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


text:
Find pattern that is present twice and allow <=2 mismatches on each
I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with <= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-

TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA

Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.

Do you have any additional assumptions you can make? For example, are the pairs of barcodes within a specific distance from each other  or the reads ends? I ask because with an edit distance of two (including indels) the chances of random matches will start to increase

This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.
I just copy pasted your examples seqs into files to read (called them seqs and code).
#need blast installed
import subprocess
import pandas as pd

path = '/Users/ur/dir/'
#make a blast db
with open (path + 'seqs.txt', 'r') as fin:
    with open (path + 'seqs.fa', 'w') as fout:
        for i, line in enumerate(fin):
            fout.write('>'+str(i)+'\n')
            fout.write(line)

subprocess.run([
    'makeblastdb',
    '-in', path + 'seqs.fa',
    '-dbtype', 'nucl'
])

#use blast
subprocess.run([
    'blastn',
    '-task', 'blastn-short',
    '-query', path + 'code.txt',
    '-db', path + 'seqs.fa',
    '-out', path + 'out.txt',
    '-outfmt', '6'
])

#grab seqs with 2 hits (ie ~2x len(query)
df = pd.read_csv(
    path + 'out.txt',
    sep='\t',
    header=None,
    names=['query', 'db_name', 'percent','len'],
    usecols = [0,1,2,3]
)
df = df.groupby('db_name').agg('sum')
print(*list(df.query('len > 27').index))

I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.


text:
Comparison of gene set enrichment statistics
I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.
As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.
How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply.
I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:
library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")

# Create a design matrix and
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)),
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
##
## $correlation
## [1] -0.02036983
##
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
##
## $correlation
## [1] 0.3024643
##
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
##
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance.
As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:
               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value?

I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications.

Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used inter.gene.cor = NA in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower.
Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (cameraPR) because it doesn't have the same results as in camera and it doesn't allow to estimate the intercorrelation.

Noticed an error in your code - in the call to topTable, you need to include 'sort.by = "none"' or else it reorders by significance, making the 1:40 and 41:80 indices for geneSetTest no longer applicable.The geneSetTest p-values for the three gene sets now come out to 1, 9.04e-13, and 0.97. It seems that the geneSetTest p-value is arguably more "correct" here, since only GS2 truly shows differential expression between samples 1-5 and 6-10. CAMERA says they're all significant!

@NuclearWang Nice catch, that's why I usually use ids2indeces! I am not sure if this is the real case, I'll try to expand the answer using the ALL dataset in a couple of days. Hope this helps! If you find a better answer, please let me know, I am interested in those methods :)

Why is CAMERA giving a significant p-value for the first comparison?

@flies There is only one comparison. Do you refer to the enrichment in the GS2?

@Llopis I refer to GS1. AFAICT, rows 1 through 40 are just a uniform normal distribution across samples, so the test should be negative.

@flies Yes, you are right. The GSEABenchmarkeR aims to find a consensus. My understanding is that the low correlation and the VIF makes strange things. But it might be an effect of a small dataset...

@Llopis thanks for your response on this old question! It's really odd and rather disturbing that this is what you found. I will say that I recently had a case where I used cameraPR where the ids2indices was scrambled with respect to the expression input (I used t-statistics output by limma), and the significance was all negative. Once I fixed the bug that scrambled the indices, I got good results that made biological sense. So camera will give negative results in some cases - at least with fixed inter.gene.cor!


text:
When working with Bioconductor and the curatedOvarianData dataset, what are the unit values of the expression data?
I am wondering how to interpret values inside the gene expression levels contained within Bioconductor. For example, if we have the following commands via R, we get:
source("https://bioconductor.org/biocLite.R")
biocLite("curatedOvarianData")
library(curatedOvarianData)
data(TCGA_eset)
head(exprs(TCGA_eset)["CXCL12",])

TCGA.20.0987 TCGA.23.1031 TCGA.24.0979 TCGA.23.1117 TCGA.23.1021 TCGA.04.1337
    5.437923     6.307250     4.193748     9.071293     5.101532     6.697770

The above calls out gene expression levels for the gene CXCL12. However, I am wondering how to interpret the values 5.437923, 6.307250, etc.
I know they are expression levels, but do they come with units? Do larger values mean greater cancer possibility?

The values are rma normalized expression values from affy microarrays. They are indeed log2 transformed (as @gringer already expected). You can find more info about the data here, under heading TCGA_eset.
To get more info about rma, check this out.
Larger values are not necessarily greater cancer possibility, I wish biology was that simple.

These values represent the expression of a particular gene in each sample. While expression of an immune system gene like CXLC12 might be associated with the aggressiveness of a cancer, that's not going to be the case for the majority of genes.
Looking at the numbers, these look like they are probably log$_2$-transformed values. In order to know for sure, plotting a density distribution of the dataset is a good idea:
plot(density(exprs(TCGA_eset)));

If the values are in the range of around 0..15, then it'll be log$_2$-transformed values, whereas if they're in the range of 0..100000, it'll be "raw" expression values. The log-transformed values are preferred for data analysis because gene expression tends to be close to a normal distribution in log-space.
These numbers are best treated as relative counts (i.e. the most appropriate unit is "expression"), as variation in the preparation of the sample (or the sequencer, or the microarray) can influence the results that come out the other end. Converting them to something more meaningful (e.g. "transcripts per cell") requires normalising against a known reference sample.


text:
Difference between 1D and 1D^2 data
I am working on ONT data. Initially, I have worked on data from 1D using Minimap2 aligner. I came across 1D^2 from ONT website. I was wondering how different is the data from two techniques? What is the difference between 1D and 2D? and do we get 2  fastq files (one per strand) for each sample?
TIA
I would greatly appreciate any help regarding my question.

1D² data that is called by Albacore generates two overlapping sets of files. Here's what Devon O'Rourke mentioned on the ONT community forums and on Twitter:

Following Albacore basecalls on a 1d2 library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D2 script. I'd like to use the reads generated from these scripts in a genome assembly.
After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D2 directories - namely, those 1d reads that generated 1d2 reads!

Link here, for those who have ONT community access.
I've recommended generating a list of the 1D^2 reads, then excluding those reads from the 1D fastq files. I have a couple of scripts that can help for this:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

More information on this can also be found in this StackExchange answer.


text:
getting KeyError on certain residues in BioPython: ATOM vs. HETATM lines
When I want to load atoms for residuum (' ', 85, ' ') in chain A for protein 3tmm with structure[0][chain_id][residue], but I'm getting KeyError exception:
Traceback (most recent call last):
  [...]
  File "/home/pok/Checkouts/bilkoviny/loader.py", line 118, in _get_residue_atoms
    return list(self.structure[0][chain_id][residue])
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Chain.py", line 74, in __getitem__
    return Entity.__getitem__(self, id)
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Entity.py", line 40, in __getitem__
    return self.child_dict[id]
KeyError: (' ', 85, ' ')

I assume this is because this residuum do not have "ATOM" lines in PDB file:
[...]
ATOM    704  CD1 ILE A  84      19.348  -4.839  18.138  1.00 47.08           C
ANISOU  704  CD1 ILE A  84     5862   6108   5917   -216  -1513    609       C
HETATM  705  N   MSE A  85      21.619  -0.130  16.236  1.00 31.89           N
ANISOU  705  N   MSE A  85     5143   3072   3900   -766  -1011    962       N
HETATM  706  CA  MSE A  85      22.277   0.675  15.226  1.00 40.76           C
ANISOU  706  CA  MSE A  85     6708   3869   4911   -972   -729   1162       C
HETATM  707  C   MSE A  85      23.146   1.786  15.818  1.00 44.98           C
ANISOU  707  C   MSE A  85     7191   4141   5757  -1268   -558   1109       C
HETATM  708  O   MSE A  85      24.257   2.045  15.338  1.00 43.69           O
ANISOU  708  O   MSE A  85     7055   3855   5692  -1611   -151   1153       O
HETATM  709  CB  MSE A  85      21.244   1.265  14.269  1.00 46.50           C
ANISOU  709  CB  MSE A  85     8043   4309   5316   -709   -948   1345       C
HETATM  710  CG  MSE A  85      21.849   1.954  13.052  1.00 54.03           C
ANISOU  710  CG  MSE A  85     9604   4929   5995   -907   -623   1638       C
HETATM  711 SE   MSE A  85      20.503   2.755  11.877  1.00 84.35          SE
ANISOU  711 SE   MSE A  85    14348   8336   9366   -486  -1038   1884      SE
HETATM  712  CE  MSE A  85      21.690   3.538  10.518  1.00 83.79           C
ANISOU  712  CE  MSE A  85    15063   7865   8907   -915   -414   2326       C
ATOM    713  N   SER A  86      22.636   2.453  16.846  1.00 43.76           N
ANISOU  713  N   SER A  86     5234   6599   4794  -1778   -946    195       N
[...]

I have been told "MSE is same as MET, it is just with selenium instead of sulphur" and that I probably need to load these atoms (I'm approximating secondary structure element position, so I want individual atom coordinates).
Is there some way how to get atoms for this residue in this protein via BioPython?

You can find the answer in The Biopython Structural Bioinformatics FAQ:

What is a residue id?
This is a bit more complicated, due to the clumsy PDB format. A
residue id is a tuple with three elements:

The hetero-flag: this is 'H_' plus the name of the hetero-residue (e.g. 'H_GLC' in the case of a glucose molecule), or 'W' in the case
of a water molecule.
The sequence identifier in the chain, e.g. 100
The insertion code, e.g. 'A'. The insertion code is sometimes
used to preserve a certain desirable residue numbering scheme. A Ser 80
insertion mutant (inserted e.g. between a Thr 80 and an Asn 81
residue) could e.g. have sequence identifiers and insertion codes as
follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering
scheme stays in tune with that of the wild type structure.

The id of the above glucose residue would thus be ('H_GLC', 100, 'A').

So in your case it would be ('H_MSE', 85, ' ').

Thank you! Please do you know if there is a way how to get something like (*, 85, *) - I mean, if I do not care about hetero-flag and insertion code at all.

@jhutar ignoring biopython hetero-flag would make sense, as it is biopython-specific workaround for incorrect PDB files; I don't know if it's possible. Do not ignore insertion codes, they are part of the sequence number.


text:
how to add a custom elements into PyMOL visualization?
I would like to add some 3D elements (balls and lines) to the PyMOL visualization to highlight some "holes" in the structure. Say I fetch 1tqn and I want to add a yellow ball somewhere (to specific coordinates, say [-19, -23, -14]) in the view. Or a bar from - to coordinates.

I was going through wiki, did a online search but have not found a way how to add these simple 3D entities into the scene. I guess it should be trivial and I'm just using wrong queries. Because I need to add lots of these, I need to be able to do that via scripting interface.
Thank you in advance for any help!

see here https://stackoverflow.com/questions/2060582/draw-a-colored-sphere-from-cartesian-coordinates-in-pymol/76749167#76749167 Draw a colored sphere from cartesian coordinates in pymol  and here https://pymolwiki.org/index.php/Load_CGO and https://pymolwiki.org/index.php/CGO_Shapes

my attempt :
import pymol

from pymol.cgo import *

from pymol import (

                    cmd ,
                    )

print('########## PYMOL VERSION ##########################################')
print('         ',  cmd.get_version() )
print('###################################################################')

pymol.finish_launching()

cmd.fetch('1tqn')

cmd.remove('solvent')

pos = cmd.get_position()

r = 3

print('###############################################################################\n')

print('### POSITION : ' , pos,'  ###\n')

print('###############################################################################\n')

spherelist = [
    [COLOR,    1.000,    1.000,    0.000,
    SPHERE,   pos[0] , pos[1] , pos[2]  , r ],

    ]

for i in spherelist:

    a = str(spherelist.index(i))

    cmd.load_cgo(i, 'sphere_'+a )

cmd.zoom()

cmd.set('cartoon_transparency' , '0.5')

result:

Aaah, load_cgo() https://pymolwiki.org/index.php/Load_CGO ("Compiled Graphics Objects"). There are also other shapes: https://pymolwiki.org/index.php/CGO_Shapes. Thank you!


text:
Empty .result in MirDeep Star
I am using MirDeep star tool for RNA sequence alignment. I have tried both GUI version and command line version of this tool. But every time, I am getting empty . result file and .known_mir file. This tool is neither giving any error nor any warning with my data. I have also tried the demo file given with this and this tool is working fine with this (i.e. giving some results in .result and known_mir file).
Can anyone suggest me what could be the possible reason for this. Why I am getting empty file.

I have found reason for not getting any result. I was using default adapter given by mirdeep* which gives me zero alignment with my data. This error pointed out by author of mirdeep* ( thanks to Jiyuan and colleen). After giving right adapter, I got the right result.


text:
Filter reads belonging to unique sequences with threshold
I have CLL samples as fastq files and I want to remove those reads which have a unique sequence with less than 10 read counts each. I tried this by following way using awk to make it faster:
for fq in $(ls fastq1/*.fq.gz); do
    basename=$(basename "$fq" | cut -f1 -d '.')
    zcat $fq | awk "/^@/{getline; print}" |sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt; grep -B 1 -A 2 -f tmp.txt $fq | grep -v ^--'| pigz -p 10 > fastq2/$basename.fastq.gz

done

rm tmp.txt

But this is not working well as many reads with more than 10 read counts are also missing. I think it's OK to stick with bash/awk for fast fastq file processing. What am I missing here, or is this command wrong. Can anyone suggest me the right and fast way to do it??

I don't know fastq, but at the start of the pipeline you zcat $fq and in the middle you grep -B 1 -A 2 -f tmp.txt $fq -- do you need to zcat $fq | grep -B1 -A1 -f tmp.txt ?

Either zcat | grep when you're trying to pick reads based on sequences, or zgrep. This code should not work as-is.

@glennjackman: In tmp.txt, I am collecting all unique sequences and finally using this file to count the number of occurrences in the original file. I tried to make it continuous without creating any tmp.txt but that didn't work. Rest other arguments of grep are used to take a complete chunk of 4 lines in fastq files it satisties the criteria of having the number of occurrences more than 10. For more information, you can check the fastq file format.

@RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python.

Your file $fq is either gzipped or not gzipped. If it is gzipped, grep ... $fq should not work. If it's not gzipped, zcat $fq should not work. I don't think both of those statements can work with the same $fq file.

My file $fq is surely gzipped. For this reason, I started with zcat. As per the answer suggested below, zcat|grep can work well with extracted unique sequences in the second step.

You say that you are "trying awk to make it faster". What else have you tried that was slower?

I have not encountered this problem before, and I am wondering if this is an XY problem. Why are you trying to do this?

Checking for line start to find the read header is not the best option. Even if unlikely it is possible that the qualitys starts with a @. It's better to compare the current line number. Also you can use awk directly to count and print the result.
zcat $fq|awk 'NR%4==2 {count[$0]++} END { for(sequence in count) { if(count[sequence]>=10) print sequence }}'

As the other statet out, the second command should be fine if you use zcat instead of cat.

Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.
I'd adjust your script to the following:
for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
 zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

Instead of

for fq in $(ls fastq1/*.fq.gz)
use this
for fq in fastq1/*.fastq.gz
Iterating over ls is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.
Rather than removing the separators that grep will introduce by filtering them out, use the following option to not create them:
grep --no-group-separator
And, minor bug, as others have mentioned, you should use zgrep or zcat | grep since your files are compressed.

Oh - and quote your variables.
Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.

Thanks for your response. I tried zgrep, but I am not getting any output. Although the script is working till generating tmp.txt. But after that, I didn't see any output file in fastq2. So I tried zcat|grep and it was working. I'll check whether it is giving a correct response or not after completing the pipeline.

@Lot_to_learn  zcat | grep may be a better option depending on your setup, as the --no-group-separator option may not work properly with zgrep


text:
Somatic mutations for normal WES samples
I am trying to identify somatic mutations on healthy subjects (specifically in WES samples). As of now, I have found a couple of papers who have done similar studies such as this paper where somatic mutations in healthy subjects are identified (using RNA-SEQ data). But in my case, I want to find somatic mutations from normal WES samples. For somatic mutations, we need both tumors and matched normal samples to identify somatic mutations but how to identify somatic mutations from normal samples only. Please suggest. Thanks.

You can follow the GATK pipeline

@csgroen: As I have mentioned in my question, I have to find somatic mutations in normal samples only. The link you have shared contains the GATK pipeline for somatic mutations with tumor and matched normal samples. I am already aware of this. Thanks.

Hm, I see... The issue is that you need your 'baseline' germline genome to be able to call somatic mutations. Otherwise, how can you know if they are germline or somatic?

Yes. Precisely. Or I just found one way in this paper (Whole-genome DNA sequencing provides an atlas of somatic mutagenesis in healthy
human cells and identifies a tumor-prone cell type) in which conditions are not clear to me.

From a quick read, my interpretation is that they aligned against a reference genome (hg19), then removed annotated known germline SNPs (using dbSNP as reference)

check out this one:https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-, this is for germline mutations

@PhoenixMu he wants somatic, not germline mutations. That's actually the challenge, as it is normal tissue and doesn't have a proper reference. But personally, I think what I described from my interpretation of the paper you shared seems to be at least on the correct path. Align against reference genome and remove common SNPs. It's not perfect, there will be false postiives and negatives, but it should work overall.

@csgroen: Yes. I agree with you. But my problem is the result's accuracy. The pipeline used in the above-mentioned paper considered 3 variant callers and allowing only variants with low allele fraction. How low value of allele fraction gives a better solution?? And 3rd condition "variants that occurred during the in vitro culture of
the clone" is also unclear to me.

Paper quote: "To eliminate all other germline variants, including rare ones, we utilized the low confidence germline variants called by GTEx. [...] Our goal was to exclude as many germline variants as possible, so we reasoned that using all germline variants called by GTEx—including the low-quality ones—was the safest option to minimize germline mutation contamination in our somatic mutation calls." I believe they removed all annotated variants, even if they were low confidence. I don't understand your in vitro comment though, makes no sense for this study.


text:
Add tags to the read pairs in a bam file
I have a bam file from paired-end sequencing, in which each R1 has the RX tag assigned (via UMItools group) denoting a Unique Molecular Identifier (UMI).
I want to call molecular consensus reads using fgbio but fgbio GroupReadsByUmi gives me an error, because it expects both R1 and R2 to have the RX tag set.
I would like a tool that will go through my bam file and assign each R2 the same RX tag as the R1.

For something that specific you'll probably need to write something.

I would also ask this of the umitools developers: https://github.com/CGATOxford/UMI-tools/issues

The best thing would be to figure out how to do this in UMItools. I don't know that tool at all, but it seems that it is supposed to be able to deal with / tag both reads in paired data. e.g. did you use --paired flag, etc.
I would suggest investigating that route first before you look at the next part of my answer.

Per Devon Ryan's comment suggestion I would suggest a workflow similar to this:

sort bam file by query name (samtools sort -n) to get the R1s immediately following the R2s
Write a short script/command that

reads through the bam file and for each line
looks up the tag in the R1, prints R1 line
reads the R2, substitutes the RX tag in, prints R2 line

pipe that through samtools view and possibly samtools sort back to a bam file of the same type that you started with.

My pysam is a little rusty but something like this might work
import pysam

bam = pysam.AlignmentFile("bam", "rb")
# print header
print(bam.header)

tag = None
for read in bam:
    r2 = False
    if tag:
        read.assign_umi_tag() # assign to read- not sure where/how
        r2 = True
    else:
        tag = read.lookup_umi_tag()  # i don't know where this tag is in bam format
    print(read)
    if r2:
        tag = None

I haven't tested that in any way, caveat emptor.
I am not personally familiar with UMItools and how it assigns tags so can't suggest anything more specific, sorry.

Thanks. My pysam is non-existent, but I can probably do it in Rsamtools.


text:
How to anonymise a bam file
I have a bam file of, for instance, RNA-Seq, which contains patient-identifiable data in the form of Single Nucleotide Polymorphisms (SNPs) throughout.
I would like a technique to take this aligned bam file and replace all the non-reference bases with the reference bases at the equivalent positions (hard-clipping any soft-clipped bases while doing so), in order to have an anonymised bam file without the patient-identifiable information.
What is the best way to do this? I want the resulting anonymised bam file to have the reads in the same positions, for downstream tools that care about the alignment positions.

Hi @KraZug normally patient confidentiality is maintained through anonminity of name and personal information, but not to anonymize the DNA by ad hoc editting. No ethics committee would request that. Personally, I don't agree here.

In the event of a publication, normally you can submit to ega or dbgap https://www.ncbi.nlm.nih.gov/gap/docs/submissionguide/#:~:text=To%20comply%20with%20HIPAA%2C%20personally,submitted%20in%20dbGaP%20data%20files. and they have a set of guidelines on how u de-identify samples

You can simply replace the sequence and quality scores with a single * in each field. The CIGAR information will still be there, so it will be apparent where there were InDels and possibly mismatches (assuming your aligning uses = and X CIGAR operations), but since the actual sequence of these variants cannot be known this is quite anonymized.


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:
Sort vcf by contig and position within contig
Due to tabix constraints, I need to sort a vcf so that contigs and then positions within contigs occur in numerical order in the vcf.  I don't know if the following will sort positions within contigs, but I tried this:
cat in.vcf | vcf-sort -c > out.vcf

but was informed:
Old version of sort command installed, please run without the -c option.

I've read that can cause problems here. Further, my contigs are't numeric, but are something like this: 'Erumb1_s00000030'.
I see two possible solutions.  (1) I am on Mac OS Mojave, and if there were a way to install the a new unix-style sort, I could use it. I tried homebrew and didn't find anything. (2) there are likely some bash commands that would do it.  If anyone has an idea in either of these directions, I would appreciate it.

It might be worth posting a question about these tabix constraints. We could be able to help you fix that which would allow you to use bcftools.

You could simply use bcftools sort for it:
$ bcftools sort input.vcf > output.vcf

If you really want to use bash only, you can do this:
$ grep "^#" input.vcf > output.vcf
$ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf

The first command will write the header information to the new vcf file.
The second will sort by contig name and position and append the result to the new vcf file.
fin swimmer

I get this: Writing to /tmp/bcftools-sort.7BfymN
[W::vcf_parse] Contig 'Erumb1_s00000030' is not defined in the header. (Quick workaround: index the file with tabix.).  So, this is where I started, having to index the .vcf.  Looks like the grep approach worked, but I've got to give it a closer look.

@PeterPearman first check that your sort implementation deals with the -k flag correctly. The command in this answer should work fine for GNU sort (the default on Linux) but I'm not sure it works correctly with BSD sort (default on macOS). The man page I found seems to suggest that the V and g modifiers aren't supported. However, you don't really need the V and you can use n instead of g: sort -k1,1 -k2,2n.

bcftools is very strict about the correct vcf structure. It heavily relies on the header data. This should include contig information and any key used in the INFO and FORMATcolumn. You should always have this information there. bcftools is a very powerful tool. So be strict to your vcf file, that you can use it.

@terdon is right. V isn't needed. I'm used to it, handling human chromosome names. Doing it that way chr2 is located before chr10. But in vcf file the only conditions are, that the contig names build blocks and inside this blocks the position is ascending.

Use a subshell. First print the header and then the variant lines, sorting by chromosome and position:
(grep ^"#" input.vcf; grep -v ^"#" input.vcf | sort -k1,1 -k2,2n) > sorted.vcf


text:
SNPEFF genome build with specific codon usage
Failing to build my own Arabidopsis thaliana reference with Mt using 'codon.Mitochondrial' and Pt using 'codon.Bacterial_and_Plant_Plastid' but the same commands works well with both lines with codon tables removed below.
The snpeff commands
java -Xms24g -Xmx24g -jar /opt/biotools/snpEff/snpEff.jar build -gff3 -v Arabidopsis_thaliana.TAIR10.46

returns:
java.lang.RuntimeException: Error parsing property 'Arabidopsis_thaliana.TAIR10.46.Pt.codonTable'. No such codon table 'codon.Bacterial_and_Plant_Plastid'
        at org.snpeff.snpEffect.Config.createCodonTables(Config.java:169)
        at org.snpeff.snpEffect.Config.readConfig(Config.java:650)
        at org.snpeff.snpEffect.Config.init(Config.java:480)
        at org.snpeff.snpEffect.Config.<init>(Config.java:117)
        at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
        at org.snpeff.SnpEff.run(SnpEff.java:1183)
        at org.snpeff.SnpEff.main(SnpEff.java:162)

(same error for the Mt when removing Pt)
my fasta has the two chromosomes
>1 dna:chromosome chromosome:TAIR10:1:1:30427671:1 REF
>2 dna:chromosome chromosome:TAIR10:2:1:19698289:1 REF
>3 dna:chromosome chromosome:TAIR10:3:1:23459830:1 REF
>4 dna:chromosome chromosome:TAIR10:4:1:18585056:1 REF
>5 dna:chromosome chromosome:TAIR10:5:1:26975502:1 REF
>Mt dna:chromosome chromosome:TAIR10:Mt:1:366924:1 REF
>Pt dna:chromosome chromosome:TAIR10:Pt:1:154478:1 REF

my data folder is OK
data/Arabidopsis_thaliana.TAIR10.46/
drwxr-xr-x   3 u0002316 domain users 4.0K Apr 15 13:33 .
drwxr-xr-x 244 u0002316 domain users  16K Apr 10 09:39 ..
-rw-r--r--   1 u0002316 domain users 106M Apr 10 09:40 genes.gff
-rw-r--r--   1 u0002316 domain users 117M Apr 10 09:40 sequences.fa

my snpeff.config contains:
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : codon.Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : codon.Bacterial_and_Plant_Plastid

the two codon tables are present by default in the snpeff.config
codon.Mitochondrial                         : TTT/F, TTC/F, TTA/L+, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/W, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

codon.Bacterial_and_Plant_Plastid           : TTT/F, TTC/F, TTA/L, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/*, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

Anybody ever succeeded to build a custom genome with specific codon tables?
Thanks in advance

codon table issue
I solved the issue with codon tables, the correct way to specify a table is the following:
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : Bacterial_and_Plant_Plastid

Note that one should remove the leading prefix 'codon.' from the name of the tables chosen from the config file
Concerning accessory files:
data to reproduce this issue:
baseurl="ftp://ftp.ensemblgenomes.org/pub/plants/release-46"
wget ${baseurl}/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ${baseurl}/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gff3.gz
wget ${baseurl}/variation/vcf/arabidopsis_thaliana/arabidopsis_thaliana_incl_consequences.vcf.gz

I suspect that the format used by Ensembl is not 100% OK for SnpEff; I therefore simplified the fasta headers to keep only the first argument (name) and remove the remaining of the line
eg.
>AT3G05780.1 pep chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 transcript:AT3G05780.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]

becomes
>AT3G05780.1

Despite this modification, I still have an issue with
java -jar /opt/biotools/snpEff/SnpEff.jar build -gff3 -v Arabidopsis_thaliana_TAIR10_46 2>&1 | tee -a log.txt

a LOT of warning about not being able to add UTR while it seems that UTR's are well present in the GFF3 file

Reading GFF3 data file  :
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  WARNING: Cannot find transcript 'transcript:AT1G01010.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Cannot find transcript 'transcript:AT1G01020.3'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Cannot find transcript 'transcript:AT1G01020.5'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Cannot find transcript 'transcript:AT1G01020.4'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Cannot find transcript 'transcript:AT1G01020.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Cannot find transcript 'transcript:AT1G01030.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line 155  '1  araport11   three_prime_UTR 11649   11863   .   -   .   Parent

nothing good coming out of the cds.fa input although I treated it with the same name fix as the pep file

00:00:19  CDS check (FASTA file):
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'
00:00:19  Checking database using CDS sequences 00:00:19  Reading CDSs
  from file
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'...
  00:00:20  done (75976 CDSs). 00:00:20 Comparing CDS...    Labels:         '+' :
  OK        '.' : Missing       '*' : Error
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
many more lines with dots
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
    ...............................     CDS
  check:    Arabidopsis_thaliana_TAIR10_46  OK: 0   Warnings: 0 Not found:
  104147    Errors: 0   Error percentage: NaN% Fatal error: No CDS checked.
  This is might be caused by differences in FASTA file transcript IDs
  respect to database's transcript's IDs. Transcript IDs from database
  (sample):     'transcript:AT3G12250.7'    'transcript:AT5G18570.1'
    'AT3G27810.1'   'AT3G27810.2'   'transcript:AT2G16810.1'
    'transcript:AT3G01170.1'    'AT4G37030.1'   'transcript:AT1G79940.4'
    'transcript:AT2G44850.1'    'AT1G49540.1'   'AT1G49540.2'
    'transcript:AT4G35290.5'    'transcript:AT4G15840.2'
    'transcript:AT4G27500.1'    'AT4G37022.1'   'transcript:AT2G40008.1'
    'AT3G27820.1'   'transcript:AT4G24265.1'    'transcript:AT2G20618.1'
    'transcript:AT4G18375.2'    'AT4G37020.2'   'AT4G37020.1' Transcript IDs
  from database (fasta file):   'AT4G02541'     'AT1G12410.1'   'AT4G26510'
    'AT4G02540'     'AT3G04500.1'   'AT3G04500.2'   'AT3G04500.3'   'AT2G36250'
    'AT5G63350'     'AT2G36255'     'AT2G12280'     'AT1G03970.1'   'AT5G18390.1'
    'AT3G04490.4'   'AT5G38386'     'AT3G04490.2'   'AT3G04490.3'
    'AT3G04490.1'   'AT5G38380'     'AT4G02530'     'AT4G26520'     'AT1G69280.2'

link to the full stdout log after running the command without cds.fa
https://www.dropbox.com/s/gxa4e533p64wxm9/SnpEff_build-log.txt?dl=1


text:
Filter bam with sambamba
I have a very large bam file and I want to filter it to keep only a handfull of positions as defined in a bed file. Can I do this with sambamba or do i need another tool?
Thanks

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

I think it is possible:
sambamba view --regions regions.bed --with-header --format=bam reads.bam  > filtered.bam

Hi @Pablo could you kindly edit your question inserting this as a suggestion?

You need bedtools, specifically bedtools intersect command, not sambamba:
bedtools intersect -abam A.bam -b B.bed > A_in_N.bam


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:
R package grur: missing_visualization() fails during redundancy analysis of .vcf
I am trying to use the R package grur to process GBS SNP data, which is installed from github as instructed.  For starters, I am running the test code in this vignette.  As suggested I am running it from console in RStudio (v. 1.1.383), and with a recent R installation (below).  I do, and get:
 > ibm <- grur::missing_visualization(data = "example_vcf2dadi_ferchaud_2015.vcf", strata = "strata.stickleback.tsv")
#######################################################################
#################### grur::missing_visualization ######################
#######################################################################
Folder created:
missing_visualization_20181209@1523

Importing data

Reading VCF...
Generated a filters parameters file: filters_parameters_20181209@1523.tsv

Number of SNPs: 31802
Number of samples: 177

conversion timing: 1 sec

VCF: biallelic SNPs
Cleaning VCF sample names

Synchronizing sample IDs in VCF and strata...
Reads assembly: reference-assisted
Filters parameters file: updated

Number of chromosome/contig/scaffold: 196
Number of locus: 17095
Number of markers: 31802
Number of individuals: 177

Working time: 1 sec

Deprecated function, update your code to use: filter_monomorphic

Scanning for monomorphic markers...
Number of markers before/blacklisted/after: 31802/0/31802

Tidy genomic data:
    Number of markers: 31802
    Number of chromosome/contig/scaffold: 196
    Number of individuals: 177
    Number of populations: 8

Informations:
Number of populations: 8
Number of individuals: 177
Number of ind/pop:
HAD = 21
HAL = 19
KIB = 17
KRO = 20
MOS = 20
MAR = 20
NOR = 20
ODD = 40

Number of duplicate id: 0
Number of chrom/scaffolds: 196
Number of locus: 17095
Number of SNPs: 31802

Proportion of missing genotypes (overall): 0.022169

Identity-by-missingness (IBM) analysis using
    Principal Coordinate Analysis (PCoA)...
Generating Identity by missingness plot

Redundancy analysis...

Error in UseMethod("ungroup") :
  no applicable method for 'ungroup' applied to an object of class "character"
In addition: Warning message:
Trying to compute distinct() for variables not found in the data:
- `strata.select`
This is an error, but only a warning is raised for compatibility reasons.
The operation will return the input unchanged.

The final object is not generated, however some plots of the PCoA analysis are. I get this result on two Macs, one with Sierra, the other with High Sierra, but both with the same RStudio and R versions. All R packages are up to date. I get the same error when I run my own .vcf
Does anyone else get these errors on this test dataset from the repo? Can anyone suggest how to deal with this?
> version
platform       x86_64-apple-darwin15.6.0
arch           x86_64
os             darwin15.6.0
system         x86_64, darwin15.6.0
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

This seems to be an internal error (using dplyr::ungroup without doing dplyr::group) . If you provide a minimal reproducible example, the maintainer of the package might be able to solve it (The maintainer should be able to reproduce the error, so provide enough data to make it possible).


text:
Parsing SAM/BAM files for additional information
I used BWA-MEM to alignment and I would like to gather the some informations like total % of match, mismatch, insert/delete etc. I am wondering if there is any existing tools that produces this information?

Other options... samtools stats maybe gives you some of the stats you are after. Chances are you already have samtools so nothing additional to install. Also useful may be pysamstats, in particular the variation sub-command.

pysamstats looks really neat!

using BioAlcidaeJdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html , you can loop over the cigar strings of each mapped read and extract the number of inserted bases, clipped bases, etc...
java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar()!=null).flatMap(R->R.getCigar().getCigarElements().stream()).collect(Collectors.toMap(C->C.getOperator(),C->(long)C.getLength(),(L1,L2)->L1+L2)).forEach((O,N)->println(O.name()+"\t"+N));' toy.bam
I   15
H   11
S   1
D   1
N   14
P   2
M   191

to get the number of mismatches, you can preprocess the bam with http://lindenb.github.io/jvarkit/SamFixCigar.html

I've written my own script, readstomper, which gives per-site statistics from mpileup results. Here's an example usage:
$ samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.csv

Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00

It defaults to reporting variants as proportions of the total coverage, but this can be changed to reporting the individual counts:
$ samtools mpileup -B -Q 0 -f mmus_chrM.fa mm2_1.2M_vs_Mmus_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl --counts

Assembly,Position,Coverage,ref,skip,cR,pR,A,C,G,T,d,i,InsMode
chrM,1528,491,C,0,409,0.87,17,0,2,15,25,65,A;52
chrM,1529,492,A,0,366,0.75,0,2,2,9,112,12,A;5
chrM,1530,493,A,0,409,0.83,0,2,1,13,68,7
chrM,1531,495,A,0,463,0.94,0,1,10,0,21,8,G;5
chrM,1532,496,A,0,463,0.94,0,2,5,8,16,4,G;3
chrM,1533,495,G,0,459,0.93,3,1,0,2,30,4
chrM,1534,496,A,0,451,0.91,0,0,4,1,39,23,G;21
chrM,1535,503,G,0,391,0.78,4,1,0,4,103,9,A;6
chrM,1536,506,G,0,496,0.98,6,0,0,0,4,10,A;3
chrM,1537,506,G,0,488,0.96,8,1,0,4,5,5,A;3
chrM,1538,506,A,0,483,0.95,0,4,3,4,12,12
chrM,1539,506,C,0,466,0.92,9,0,2,14,14,27,A;6
chrM,1540,506,A,0,469,0.93,0,0,17,6,14,11,T;3
chrM,1541,506,G,0,453,0.90,12,1,0,5,31,21,A;4
chrM,1542,505,C,0,421,0.83,6,0,2,14,62,8,CT;3


text:
Filtering out 8-oxoG-derived artefacts in variant calls from cfDNA
I have cfDNA samples which were sequenced (at about 300-1000X) on a 30-gene panel via a capture-seq protocol for the purposes of calling variants on circulating tumour DNA.
However, the resulting variant calls (from both varscan and Mutect2) are dominated by transversions from C>A and G>T at low allele frequencies, which are characteristic of DNA oxidation of guanine to 8-oxoguanine, and which I suspect are artefacts.
I have matched baseline tumour samples from the same patients, sequenced in the same way, (and there are a few C>A/G>T transversions at low AFs in these samples that I can filter out by a direct allele fraction filter).
My question is this, what is the best way of filtering out the false positive mutations derived from the 8-oxoG formation and retaining the true positives (which will be at low allelic fraction)?
Some background: I do have serial (cfDNA) time points for some patients, so I am loath to just look at mutations called in the baseline tumour samples, as any (true) novel variants appearing later in the time points would be missed.


text:
Customizing colors in geom_arc_bar (ggforce)
I am trying to plot scatter plot of pie charts using ggplot2 and ggforce.
Following is my code
library(data.table)
library(ggplot2)
library(ggforce)
library(dplyr)

df <- fread("data_frame_pie_plot.csv", header = TRUE,
        data.table = FALSE, stringsAsFactors = FALSE)

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

plot_obtained
Link_to_data Data_frame.
Can I fill the colours manually (discrete(based on discrete_fill column)
and continous (based on conti_fill column))?
I have tried the following codes.
##================================================================================##
# Code 1
my_colors <- c("red","white","blue")
ggplot(df, fill = sample_mut_colour) + req_data +
  theme_no_axes() + scale_color_manual(values=my_colors)

##================================================================================##
# Code 2
my_colors_1 <- c("red", "white")
my_colors_2 <- c("blue", "white")

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill,
               color = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 3
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 4
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill),
           palette = case_when(discrete_fill == 1 ~ "red",
                               discrete_fill == 2 ~ "black",
                               TRUE ~ "white"), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##=================================================================================##

But they do not seem to work. Any idea what is wrong in the code?

For a discrete scale your Code 1 example was nearly there.  Just use scale_fill_manual instead of scale_color_manual.  For more control over how the discrete values map to color values you can name the vector, too.
my_colors <- c("red","white","blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)
# or, if you want consistent output independent of the order of
# df$discrete_fill:
my_colors <- c("1" = "red", "2" = "white", "dummy" = "blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)

You can also just use a single call to geom_arc_bar for the whole data frame (rather than a list of call results from each row) and you'll get the same effect.
ggplot(df) +
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie") +
    theme_no_axes() + scale_fill_manual(values=my_colors)


text:
Picking number of samples for clustering
I have expression data for 100 samples and clinical data for 50 among them. I want to carryout k-means clustering using expression data and then do a survival analysis (Kaplan-Meier analysis) using the labels obtained from k-means.
Can I used the expression data of all 100 samples to carryout clustering or should I use only 50 samples for clustering?

Hi @user98059 I understand your question perfectly at a analytical level, but I think it would be preferable to edit it and provide a biological explanation regarding survival analysis.I'm not sure it will attract many responses. ... anyway. In summary, I would go for 100. It also depends on whether K-means is inside or outside ML methods.

In case there is no systematic difference between the samples with or without the clinical data (i.e. one set is microarray and the other is RNA-seq), I would use all of the samples. And if the cluster labels obtained with "clustering 50 samples" and those obtained with "clustering 100 samples" does not differ too much, then doing either of two would not have much effect anyway. If you would like to compare cluster labels in that way and assuming you would want to test a number of ks, it would be better to compare each combination of ks and number of samples.

Hi @Haci The thing that bothers me is it isn't clear with K-means is within supervised or unsupervised learning. In my experience it can be used for either. These are quite different calculations. At a guess it is unsupervised because in supervised learning all the samples would be used without any question.


text:
Grouping x-axis in scatter plot
I am trying to plot one scatter plot for each group, where each group has continous values in x-axis.
Reference Data
Gives continuous x-axis but all groups are combined.
ggplot(data_for_plotting,
aes(x = log10(Count), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

Gives grouped x-axis but the variation in x-axis is lost.
ggplot(data_for_plotting,
aes(x = factor(Clusters), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

Does achieve the grouping but the gives multiple plots.
ggplot(data_for_plotting,
aes(x = log10(Count),y = log10(Nonsilent.Mutation.Rate),color=factor(Clusters))) +
facet_wrap(~ factor(Clusters), shrink = TRUE) +
geom_point(alpha = 1, size = 1) + theme_bw()

I want x-axis to be grouped based on Cluster where each cluster has continous x-axis values.

You are asking for repeating values in the x-axis. That's the only way you'll maintain a continuous x axis while seeing separate lines for each Cluster. The facet is the right solution - they're not separate plots, they're multiple facets of the same plot.

I'm not sure what you're asking for. The first and the third both seem acceptable from a generic visualization point of view. "Gives multiple plots" seems to be what you want based on the prompt? Ram seems to understand you, maybe it's just me.

@MaximilianPress it's not just you. I'm also having trouble understanding the question.

Possibly draw a cartoon of what you would like to accomplish?

You have two different types of values you want to display on the x-axis: Cluster (a categorical variable) and Count (a continuous variable).
As Ram RS suggested, the easiest way to accomplish this is by taking your first plot (which uses a continuous x-axis scale) and facet along your discrete variable.
This gets you what you want, and the rest of the work is adjusting the theme to make it look how you want.
I'm guessing based on your data that you're trying to create something that looks like this figure from Alexandrov et al., Nature, 2013:

You can accomplish this by customizing the facet and the theme.
Here's an example.
gg <- (
  ggplot(
    data = data_for_plotting,
    mapping = aes(
      x = log10(Count),
      y = log10(Nonsilent.Mutation.Rate),
      color = factor(Clusters)
    )
  )
  + geom_point(alpha = 1, size = 1)
  + facet_wrap(
    ~ Cluster,
    # this ensures that all facets are shown on the same row instead of wrapping
    nrow = 1
  )
  + theme_bw()
  + theme(
    # this removes the spacing between facets
    panel.spacing = unit(0, "lines"),
    # this removes the border around facets
    panel.border = element_blank()
  )
)

There are more parts of the theme you can adjust to your liking, like the facet header (e.g. strip.text).
See ggplot2's reference documentation for details.


text:
Generating contrast matrix for limma in loop
I am trying to generate contrast matrix for different conditions using for loop.
cluster_len <- 3
req_contrasts <- c()

for (m in 1:cluster_len) {
  for (n in 1:cluster_len) {
    if (n > m) {
      req_contrasts <- c(req_contrasts, (paste0(paste0("Group", n), "-", paste0("Group", m))))

    }
  }
}

req_contrasts returns
"Group2-Group1" "Group3-Group1" "Group3-Group2"
The contrast that needs to be passed is
"Group2-Group1","Group3-Group1","Group3-Group2"
I have tried the following code to combine each string and form a single vector
# To separate each condition by a double quote and comma
temp_1 <- paste0(req_contrasts, collapse='","', sep="")
temp_1
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"
# To remove the escape character
temp_2 <- gsub("\\", "", temp_1, fixed = TRUE)
temp_2
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"

This does not work
contrast.matrix <- makeContrasts(temp_2, levels=design_mat)
Error in parse(text = ej) : <text>:1:14: unexpected string constant
1: Group2-Group1","
^

Limma throws the following error in presence of escape character. I am assuming it is because it considers that the string is yet to end.
Following input works
contrast.matrix_2 <- makeContrasts("Group2-Group1", "Group3-Group1", "Group3-Group2",                 levels=design_mat)

In my analysis the number of clusters to be analysed is 16. Hence, the contrast conditions cannot be generated manually.


text:
Declaring an array of strings and integers
I have two csv file - one is a list of strings and other has both strings and floating values. If the strings in both the files match I need to copy the complete row from the second file to a temporary array. I have declared a matrix with n rows and m columns and saved is as dataframe temp <- as.data.frame(matrix(nrow = n, ncol = m)). I am using for loop to iterate over all rows and strcmp to compare the strings. But when the data is written to temp the strings are stored as integers. How should I declare this array? Is it possible to write this loop without declaring an array in R?

This question is unrelated to bioinformatics. You should ask somewhere more appropriate (e.g. stackoverflow). You might also want to add a minimal working example, it will be easier for other people to help you (see: https://stackoverflow.com/help/mcve). Good luck !

I'm downvoting because this isn't really a bioinformatics question, but I suggest you take a look at dplyr::semi_join

The merge function of R might be useful. The function expects two data frames, with [optionally] columns to compare, and it will output all fields from both data frames where the comparison columns match.
> a <- data.frame(test=1:10);
> b <- data.frame(test=c(sample(1:10,5),sample(c("A","B","C","D","E"))),
                  output=rnorm(10));
> merge(a, b, by="test");
  test     output
1    1 -1.3800443
2    2  0.2934274
3    3  0.0868700
4    9  1.1766914
5   10  0.3793448

The merge function can be customised quite a bit (e.g. if you want NA set for non-matching rows), see ?merge for more details.

When you read in your two CSV files, you can explicitly declare the variable types for each column:
> a <- read.csv(..., colClasses=c("character"), ...)
> b <- read.csv(..., colClasses=c("character", "numeric", "numeric", ...), ...)

You could also investigate as.is or stringsAsFactors, although this can cause problems where there are duplicate strings. Using colClasses is more explicit.
Read the help documentation in ?read.csv for a more complete explanation of these operators.
You can filter with the %in% operator in match:
> strings_of_interest <- as.vector(a$V1)
> b.filtered <- b[b$V1 %in% strings_of_interest,]

If I am reading your question correctly, perhaps you can then take b.filtered and turn that into the matrix you want.
Avoid use of for loops where you can. They can slow R down considerably. A good read on this subject is in the third chapter of Patrick Burns' The R Inferno.

Is there a reason for the downvoting?

I think it is bad for a community when downvoting without any explanation. Please explain yourself.


text:
How hard is it to clean and QC gene expression microarray data?
I am a PhD student working on developing ideas for my dissertation papers. One of my planned papers will be working with some (human) gene expression data.
I have had a class on working with gene expression data, but no practical experience. From what I remember from my project in that class, you do a bunch of pre-checks and QC, then you remove batch effects, then you summarize expression across probes and normalize. I thank that part of my class project took me about three days.
But when I talk to other people about it, they say that cleaning/QCing/normalizing gene expression data is really hard, that it will take forever to get it working, and that if I do this dissertation paper I am going to end up taking extra years to graduate.
Am I missing something here? Why do the people I talk to seem to think this is so hard and time consuming -- is it really that bad? Should I avoid having a dissertation paper that uses gene expression data if I want to graduate on time?
(The data will be from Affymetrix microarrays. There are enough samples to use two 96-well plates, so batch effects will be more of a concern than they were in my class project.)

Most of the problems with normalizing data are due to a lack of proper design. How is the experiment designed? What do you plan to do if your QC fails? Will it take months to repeat the same experiment? Or it will take 1k$ more to repeat it and 2 weeks? (Who is running with the costs?) What do you expect to learn from this experiment? Aren't there already public datasets that answer your question (I don't know what do you need to do for your dissertation paper)? Why do you do NGS instead of microarrays? (Sorry for the questions, but they will give context about why your idea)

Thanks for commenting! The samples are all collected and frozen, and haven't been run yet. There is a correlation structure in the data (we sampled the same people repeatedly over time) that might complicate things. Because the samples haven't been run yet, we're hoping we can organize our plating so that all samples from the same person end up in the same column on the same plate to minimize batch effects. We are using microarrays for cost reasons. If lots of samples fail QC, the project dies, we won't have money to do it again.

If you put all the samples of an individual on the same plate you will not be able to distinguish between batch effects and individual effects!! That's the thing that's tricky when performing those analysis (Realizing the implications and later problems you will face when analyzing the data). How many samples do you have? Is there any other condition involved (origin of the sample/hospital/date/treatment/age/sex...) that could affect your plates? If it is your first time check your design with more experienced people of your lab before running anything (or spending the samples)!

More info: We'll be looking at within-individual changes over time, and between-individual differences at the same timepoint.

Very good point on the batch vs. individual effects Llopis!!! I will be discussing this with more experienced people before this gets plated/run. We do have differences in age, sex, and race. My main concern at this point is "do I really want this to be a dissertation paper, or is it going to turn into a big mess and make me a 7-year PhD student?"... so far the two more experienced people I've talked to seem to think this should NOT be a dissertation paper, but my advisor (who hasn't done gene expression before) really wants me to do it...

I was about to write an answer to your questions, but your main question is if you should do it or not, and I can't answer that. If you want perhaps you can change the question to what aspects to you need to take into account when designing a microarray experiment (which would be a very nice question with many references) which could be answered in the format of this platform. Good luck with your disseration

It sounds like you're facing an academic, rather than a technical issue. I suggest you [edit] your question here to focus only on the technical aspects. Describe the experiment you want to run in detail and ask how long you should expect the QC step to take. The question of whether to do it and whether it should be a part of your dissertation can then be asked on [academia.se].

@terdon The main question (time scale) is mostly dependent on technical aspects of the experimental design and I do not think that other communities would be able to properly discuss it.
At the same time I agree with terdon that you should [edit] your question and add all this spicy details you wrote in the comments. Also, note that your question should be answerable, which is something I doubt a bit.

I'm voting to close this question as off-topic because it's not a real technical question

I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There are many tutorials and review papers about normalization and QC of arrays. Try there and see if you can easily understand what these papers are telling you.
Read the limma user's guide, and try to follow the tutorials therein. If you understand that quickly you can try a real experiment of your own. If you don't understand what's described in there, you might want to look for another project for your thesis.


text:
Estimating computing resources needed for a GWAS?
One of my dissertation papers is going to involve a GWAS. I have never actually run a GWAS before, and do not know how to estimate the computing resources I need for it. I asked someone on my committee who said "it's all guesswork" but I don't even know how to come up with a reasonable guess!
I need to put some estimate of how much cluster time I will need (XX cores for YY hours) in order to do a dissertation grant. I've tried Googling and searching literature for how to do an estimate of this, but I'm not finding much.
Details:

GWAS, total n=350,000, 10-15% cases
Already imputed by the time I get it, will have about 96 million variants
Will be using Plink
Cluster nodes have 24 cores, and jobs can use up to 4 gigs of RAM per core. If Plink can handle being spread across nodes, I can use up to 100 cores at once.
Time is billed in core-hours, using 20 cores for one hour is equal to using one core for 20 hours
We may want to do a few smaller GWAS where we restrict by age (eg. one of all cases and controls under 25, and one of all cases and controls who are 65+)
Need to come up with an estimate of cluster time needed for each GWAS that we might run (XX cores for YY hours)

Can you please advise me on how to estimate something like this?

Is it possible for you to do a trial run on ~1/1000 to 1/10000 of your dataset?  Runtime should just scale linearly in the # of variants.

Answer from @llrs, converted from comment:
I found an article where they said that they could do GWAS in 6000 samples in 6 hours in a regular user computer. But I couldn't find any other explanation about how does it take with other kind of data. It wasn't with plink too and I think it was a bit outdated.
Answer from @christopher-chang, converted from comment:
Incidentally, you may want to use plink 2.0 rather than 1.9 for this. It automatically makes good use of all your cores; in contrast, you’d need to start ~24 copies of plink 1.9 to saturate a 24-core machine.