text:
How to simulate nanopore reads?
I have looked already here: Tools for simulating Oxford Nanopore reads
. This doesn't answer my question, because it lists a few Nanopore read simulators, but I have specific problems with each of them.
I am trying to simulate some Nanopore reads. I could use NanoSim, but that tool seems to be very difficult to use. I only have a reference genome, but you need training data as well? anyone have any experience with this? I can only find bacterial and yeast data (and I want to try this on Arabidopsis thaliana reference)
Another promising tool I found, was a paper for SNaReSim, reference here:
https://ieeexplore.ieee.org/document/8031171. However I cannot seem to find the tool itself, which is also mentioned in the other stackexchange post I mentioned above.
Does anyone have experience simulating PacBio reads with NanoSim, SNaReSim or maybe another tool? Because I am not sure about which one to use. Others like SiLiCO or ReadSim were updated 2 years ago and I believe that Nanopore technology has improved a lot in that time.
To clarify a bit on why I want to simulate data: I just want to test several tools about SV calling and simulating for a school project. All of the points mentioned in the comments I will use for the discussion, but I wanted to focus on 1 read simulator per technology and didn't know what to use, because Nanopore is difficult to simulate in that sense. That was already known to me.

Why do you need to simulate reads? There are plenty of existing public datasets for human gDNA, parasite gDNA, human RNA, human cDNA, E. coli gDNA, and I've uploaded a couple of mouse cDNA runs.

because I want to test some SV callers, and for that I need to know exactly where my SVs are, and for that I need simulated reads that contain my known SVs

@Fini, please [edit] the question to include your last comment,  and perhaps more information, the answers might change/improve with it.

Yes, you need to know where the SVs are, but simulated reads are not necessary for that. It's much better to use annotated true data (as suggested in conchoecia's answer), because it's not possible to know all the variation in DNA that can contribute to base caller error. This is a particular issue with nanopore data, because the underlying sequencer model (current change caused by DNA moving through a pore) is very complex.

I added some lines to the question to clarify. I can use this information in my discussion. Do you possibly have a paper that claims the same thing or is this just based on your own knowledge?

It is important to train an error model as NanoSim does, as we do not fully understand the error processes involved in both the nanopore sequencing process and the basecalling process. Any sort of read simulator that does not use a species-specific model is just not going to produce realistic simulated reads.
As far as datasets go, I would recommend looking at this assembly and the associated data uploaded on EMBL:

Michael, Todd P., et al. "High contiguity Arabidopsis thaliana
  genome assembly with a single nanopore flow cell." Nature
  communications 9.1 (2018): 541.

There are several data types there, as well as an assembly to train the models on.
Is it OK to train a model using reads from one individual and a reference assembly from another individual?
Generally, it is not good to map reads from one individual to the reference genome of another individual since this introduces something called reference bias. When we align reads to a reference genome generated from the same reads, or just from other reads from the same individual we remove this problem. All true reads will map to correctly assembled parts of the genome.
Why can we train a model for nanopore reads on a reference generated from ... nanopore reads? Won't we get a bad model for indels?
Yes, well, maybe. If the reference genome has been corrected to remove indels using a tool like pilon, then you can generally presume that the reference genome has very few to no indels present as a result of sequencing error from the Nanopore/PacBio reads. As a result, the model derived from aligning reads should have the correct error profile for indels.
Again, it is important that the long-read-assembly was corrected with Illumina shotgun reads from the same individual from which the long reads were derived. Otherwise the corrected reference would have a chimaera of indels and SNPs from multiple individuals.

so to get it more clear, I should use the MinION " FASTQ files (FTP)", convert this to .fasta and use this to get a training folder using ./read_analysis.py with -r as the Arabidopsis reference genome (where should i find that?). Then use "./simulator.py" and use -c to that folder and -r to get the reads from. Does it matter that the reference I want reads from is different (in size and probably sequence) from the MinION data that I use? (as long as they are similar i guess?)

Yes, you would just use the basecalled fastq/fasta data and use the reference from the same dataset to align and make the model. I have updated my question above to address using reference/data from different individuals. I'm not sure of the exact commands to run but it sounds like you have the right idea. When navigating EMBL, there are a few levels of organization. Above, I linked to the project - the highest level in the hierarchy. Clicking on the "naviagation" tab reveals the other pages for the project including the assembly: https://www.ebi.ac.uk/ena/data/view/GCA_900303355

Thanks for the clarification. So to recap (sorry I am fairly new to this kind of data). Reference Arabidopsis =  "WGS_SET_FASTA" from the assembly link.
Basecalled fastq files Arabidopsis = MinION "FASTQ files (FTP)" from the Read Files tab. Use these two to get the training folder.
Sequence to produce reads from: My own .fasta file (also Arabidopsis, but slightly different from the ones mentioned above)

Sorry, the question just updated, i have a look

nanopore fastq (FASTQ files FTP): https://www.ebi.ac.uk/ena/data/view/ERR2173373  ... Sequel fastq (FASTQ files FTP) https://www.ebi.ac.uk/ena/data/view/ERR2173371  ... reference (WGS_SET_FASTA) https://www.ebi.ac.uk/ena/data/view/GCA_900303355 . Yep - you have the right idea w/r/t how to use all of the files.

Thanks, I have a try tomorrow. I hope someone will be able to help then if something is not working after experimenting with the program first myself. But I am much further than I was before. Thanks!

So I wanted to come back to this issue, so it could help others in the future. So these files worked, although I had many troubles working with the large dataset. Just keep in mind to reduce the dataset in some way or have a computer with high memory availability (32 gig was not enough at least).

I could use NanoSim, but that tool seems to be very difficult to use. I only have a reference genome, but you need training data as well?

Karel Břinda's fork, NanoSim-H, ships with pretrained models (or "error profiles", as they're called there), eliminating the need for training data.
Minimal usage instructions to run with the default error model (trained on 1D² reads of E. Coli on an R9 flowcell):
pip install nanosim-h
nanosim-h your-reference.fasta


text:
What are some of the major questions pursued by the field bioinformatics?
I am an occasional user of tools in bioinformatics and I am fascinated by the possibility of uncovering deeper patterns in the biological realm using the help of computer algorithms. Yet I don't know much about the current state of this science. As I am considering a post-graduate degree in bioinformatics, I am interested to dig deeper into this subject and check whether the subject indeed aligns with my interest. So I start with the basic question stated in the title: What are some of the major questions pursued by the field bioinformatics?
Thanks in advance.

Welcome to the site. I think that this question is quite broad to get an answer on this site. Give it a try on www.biostars.org, however, without knowing your interests it is hard to suggest a subfield/area.

Welcome! Bioinformatics is just a term that describes an approach to doing biology. People who use ( DNA sequencing/genomics/protein informatics + programming ) to answer biological questions are after the same answers as the the people who use ( pipettes/microscopes/molecular biology + a lab notebook ).
My best suggestion is to learn about some biological concepts that are interesting to you, see how people have applied bioinformatics to that branch of biology, and think if there is something you might be interested in doing to work on that question. It is critical to have a strong understanding of the underlying biology to get to the core of new findings using bioinformatics.
A lot of the hardcore algorithm development going on today is in genome assembly with long + noisy DNA sequencing reads, aligning long+noisy reads to genomes and trying to figure out the differences between the maternal and paternal haplotypes, working on evolutionary questions with population-level genetic variation data, and figuring out how things are related to each other with phylogenetics, to name a few.
As far as a machine-learning type approach which you seem to be hinting at in your question, one very hot topic at the moment is single-cell sequencing. This is a set of technologies that allows us to get DNA, RNA, or protein sequencing from a population of single cells to look at heterogeneity in tissues, timepoints, et cetera.
Good luck! If you have some free time on your hand this website, ROSALIND, is a fun project Euler-type website where you can learn some biology while programming bioinformatics solutions.

source: Stack Exchange

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

sort: votes


text:
Does anyone have older versions of NT and 16s RNA DB from NCBI?
I am in search of past years' versions of NT and 16s databases from NCBI.  I recently contacted NCBI about this and they replied:

Unfortunately, we do not archive past versions of the blast databases. We update them regularly and simply do not have the resources to store older versions.

Any older versions would be OK, but if you happened to have 2012 and 2007 those are preferable.

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


text:
Details of ENCODE processing pipelines
I am trying to find details for how the ENCODE ChIP seq samples were processed.
Linking to a ChIPseq samples on www.encodeproject.org (e.g. https://www.encodeproject.org/experiments/ENCSR051DXE/) offers a PDF "overview" of the pipeline which says basically nothing. The page will also offer you a map of the pipeline, which basically just has headings:

Clicking on a step gives you a popup with a link to the "Chip-seq read mapping" pipeline. This just contains the same overview PDF and a link to "ChIP-seq Data Standards and Processing Pipeline". This page contains again, a vague outline of the pipeline, a link back to the page we just came from and a list of quality control thresholds.
I can't seem to find anywhere the proper details. For example, if BWA was used, was it MEM or aln? What parameters were used for the mapping? One of the steps is read filtering: what where the filters? How were they applied? There is a link to a github page for running the pipeline on DNAnexus, but I can't make head nor tail of it.

I don't really understand it but perhaps the tests.sh it might be easier to understand for you (how they call it, what they expect and so on).

They are using bwa aln option with -q 5 -l 32 -k 2.
Ref: https://github.com/ENCODE-DCC/chip-seq-pipeline/blob/541dd361f28ef7568f3743c00dfcb882c461df2c/dnanexus/input_shield/dxapp.json line 43

Updated ENCODE ChIP-seq analysis pipeline https://github.com/ENCODE-DCC/chip-seq-pipeline2 . This also features the usage of bwa aln Ref: https://github.com/ENCODE-DCC/chip-seq-pipeline2/blob/e4d53da91f7dcf16dad984b1c9ebbbf168f8ef4b/src/encode_bwa.py line: 73.

The codes are usually available at DCC's GitHub.
You can always email them if you have more questions.


text:
Using external list of PCs for clustering
I am going to use principal components (PCs) comes from calcPCA function in URD program for clustering my cells in Seurat; So
seurat@dr$pca=URD_object@pca.load

> seurat <- FindClusters(
+     object = seurat,
+     reduction.type = "pca",
+     dims.use = 1:8,
+     resolution = 1.0,
+     print.output = 0,
+     save.SNN = TRUE
+ )
Error in GetDimReduction(object = object, reduction.type = reduction.type,  :
  cell.embeddings slot doesn't exist

There is seuratToURD function but there is not URD to seurat function I guess so I could use URD object for clustering in seurat
How I can overcome this error please? I want to see when I am using PCA from URD what would be my clustering by SNN algorithm in seurat
Thanks a lot for any help

There is a way to do this, and even better--there is documentation for how to do it! No surprise coming from the Satija Lab. In the vignette they perform multidimensional scaling, but the idea is the same. cmdscale() returns the cell embeddings. SetDimReduction() is the Seurat function you are looking for. No manual editing using @ required.
The authors use add mds embeddings, but if you want to add your own pca embeddings swap out "mds" for "pca" and mds for URD_object@pca.load.
# Before running MDS, we first calculate a distance matrix between all pairs
# of cells.  Here we use a simple euclidean distance metric on all genes,
# using object@scale.data as input
d <- dist(x = t(x = pbmc@scale.data))
# Run the MDS procedure, k determines the number of dimensions
mds <- cmdscale(d = d, k = 2)
# cmdscale returns the cell embeddings, we first label the columns to ensure
# downstream consistency
colnames(x = mds) <- paste0("MDS", 1:2)
# We will now store this as a new dimensional reduction called 'mds'
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "cell.embeddings",
    new.data = mds)
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "key",
    new.data = "MDS")

See the dimensional reduction vignette for more.

Thank you, but I have already demential reduction matrix (PCs) by another program (URD) and I want to store that in Seurat object for clustering

In SetDimReduction() replace "mds" with "pca" and mds with URD_object@pca.load.

Thank you                                                                                > pbmc <- SetDimReduction(object = seuset, reduction.type = "pca", slot = "cell.embeddings",

                    new.data = URD@pca.scores)

Show Traceback
Rerun with Debug
Error in (function (cl, name, valueClass)  :
‘cell.embeddings’ is not a slot in class “data.frame”

Try str(seuset) to ensure your Seurat object is whole. You might also RunPCA() within Seurat and return to cell.embeddings to ensure the format of your external embeddings is identical. I'd also recommend changing reduction.type="pca" to something like reduction.type="pca.external" just to ensure you aren't dealing with some protected fields in the S4 object.


text:
RSVSim insertions from chr1 to chr2
edit explaining python tag
I would still rather have a solution based on the RSVSim package in R, but while waiting for someone who might have an answer I wanted to look at other solutions as well.
If a result like what metadata(sim) produces can be achieved in python as well, I can do it without RSVSim. I need especially the numbers given there, so I do NOT need "Copied BpSeqA BpSeqB_5prime BpSeqB_3prime". The numbers need to reflect the original position of insert into chr2, because I need to reflect this back to the original (unedited) chr2.
end edit explaining pyton tag
I tried using the BioConductor forums for asking this question Question on BioConductor. However I do not seem to get an answer there and after some help with other tools on this website, i thought I could also ask the question here, in the hope that someone knows something about it.
I have a reference genome, and I need to put in some structural variations. It is really easy to do in RSVSim, and used it quite a lot already.
However, when processing INS, I found out that it just copies a DNA seq from one part, and puts it somewhere else in the same DNA seq (Bascially a DUP then), or if you specify another an extra DNAseq, to cut from, It randomly cuts from either chr1 or chr2 and puts it in the other. This creates deletions in the DNA seq I want to investigate, which is not what I want.
Does anyone know how to get this kind of result?
width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT    chr1
[2]    40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC    chr2

to this:
width seq                                                                                names
[1]    32 AAAAAAAAAAAAAAAATTTTTTTTTTTTTTTT                    chr1
[2]    48 GGGGGGGGGGGGGGAAAATTTTGGGGGGCCCCCCCCCCCCCCCCCCCC    chr2

And do it 30 times, randomly in the genome? (so take a seq from chr1 and put an insertion in chr2.
I tried this and more already for days, and I just cannot figure it out.
    library(RSVSim)

    seq_random <- readDNAStringSet("/path/random.fasta", "fasta")
    seq_ref <- readDNAStringSet("/path/reference1MB.fasta", "fasta")
    genome2 = DNAStringSet(c(seq_random, seq_ref))
    names(genome2) = c("chr1","chr2")

    y <- c(30, 30, 30, 30, 30, 50, 50, 50, 50, 50, 100, 100, 100, 100, 100, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 5000, 5000, 5000, 5000, 5000)

    length_seq = width(genome2[2])
    knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr2")
    knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr1", chrB="chr2", chrA="chr1")
    knownInsertion = GRanges(seqnames="chr1", chrB="chr2")
    sim = simulateSV(output='/folder/of/output/', genome=genome2,
                     ins = 30, sizeIns=y, regionsIns=knownInsertion)

EDIT for comment section
The simple genome example:
genome = DNAStringSet(c("AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT","GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))
names(genome) = c("chr1","chr2")

I believe I already tried those parameters. For example:
knownInsertion = GRanges(IRanges(16,25), seqnames="chr1", chrB="chr2")
sim = simulateSV(output=NA, genome=genome, ins = 1, regionsIns=knownInsertion, bpSeqSize=6, random=TRUE)

results in:
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                chr1
[2]    40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1      3   12 chr1      6   15   10  FALSE  AAAAA        AAAAAA        AAAAAA

so nothing happened at all.
If I keep it simple, you get the result I mentioned before, where it is sometimes chr1->chr2 and sometimes chr2<-chr1
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE)

> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAACCCCCAAAAATTTTTTTTTTTTTTAAAATTTTTT                                                                                chr1
[2]    40 GGGGGGGGAAAAAGGGGGGGGGGGGCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1     17   21 chr1     36   40    5  FALSE AAATTT        TTTAAA        AATTTT
2 insertion_2 chr1      1    5 chr2      9   13    5  FALSE               GGGAAA        AAAGGG
3 insertion_3 chr2     33   37 chr1     12   16    5  FALSE CCCCCC        AAACCC        CCCAAA

In this case I dont even know what happened exactly, because an uneven number of Ins, should make one of the chromosomes shorter or longer?
Finally, when setting copy on (so no cut/paste)
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE, percCopiedIns=1.00)
sim
metadata(sim)

> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    50 AAAATTTTTAAAAAAAAAAAAAAAATTTTTTATTTTTTTTTTTTTTTTTT                                                                            chr1
[2]    45 GGGGGGGGGGGGGGGGGAAAAAGGGCCCCCCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1     36   40 chr1      5    9    5   TRUE               AAATTT        TTTAAA
2 insertion_2 chr1     20   24 chr1     27   31    5   TRUE               TTTATT        TTTTTT
3 insertion_3 chr1     13   17 chr2     18   22    5   TRUE               GGGAAA        AAAGGG

Again. all the examples above have random chrA and chrB. Setting these specifically, only seems possible with the GRanges and IRanges, for which I tried tons of versions already.
It should be a very simple thing according to the manual, but actually doing it doesn't seem to get the desired results, while my version of Insertions is just the normal one in my opinion? (why would you want to create extra DELS or DUPS while making insertions?).
On their own support page, it is only possible to ask a question on the forums, but they never answer it (only look at it)

What is the problem when using the simulateSV? You don't get the desired result or something else? From the manual it seems to be possible. BTW  knownInsertion should be chromosome 2 if you want all the insertions to be into the chromosome 2. But perhaps you are missing the option random = TRUE

see my updated question. I really think it is weird that this very simple and logical thing isn't working. I already looked extensively at all the available manuals out there, but there are only 2.

If someone is able to make the result I want with python, it is also fine, but then I need the exact numbers where on the chromosome the edits WERE made (so not the numbers afterwards, so that I can reflect back to the original reference genome)

Answer eventually found on the BioConductor post:
https://support.bioconductor.org/p/114354/#114585
It basically comes down to not using RSVSim, but other R packages.
Final code snippit:
genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

sampleFromGenome <- function(genome, wd) {
  chrProb <- width(genome) / sum(width(genome))
  chr <- sample(names(genome), n, TRUE, prob = chrProb)
  start <- runif(n, 1, nchar(genome[chr]) - wd)
  GRanges(IRanges(start, width = wd), seqnames = chr)
}

n <- 2
wd <- sample(5:10, n, TRUE)
insertFrom <- sampleFromGenome(genome[1], wd)
insertAt <- sampleFromGenome(genome[2], 0)

#writeTable(insertFrom, stdout()) . # replace stdout with a file name

sequence <- genome[insertFrom]
at <- split(insertAt, seqnames(insertAt))[names(genome[2])]
atSequence <- split(sequence, seqnames(insertAt))[names(genome[2])]

replaceAt(genome[2], unname(ranges(at)), atSequence)

  A DNAStringSet instance of length 1
    width seq                                                                                names
[1]    54 GGGGGGGGGGGGGGGGGGTTTTTTTTTGGCCCCCCCCCCCCTTTTTCCCCCCCC                                                                        chr2

EDIT with RSVSim solution!
Finally found a solution based on RSVSim
genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

length_seq_chr1 = width(genome[1]) - 10
length_seq_chr2 = width(genome[2]) - 10
sizes <- c(3,3,5,5)
start_chr1 <- c("list", length(sizes))
start_chr2 <- c("list", length(sizes))
end_chr1 <- c("list", length(sizes))
index = 1
for(i in sizes){
  newstartchr1 <- sample(1:length_seq_chr1,1)
  newstartchr2 <- sample(1:length_seq_chr2,1)
  start_chr1[[index]] <- newstartchr1
  start_chr2[[index]] <- newstartchr2
  end_chr1[[index]] <- newstartchr1 + (i-1)
  index <- index + 1}
start_chr1 <- as.integer(start_chr1)
end_chr1 <- as.integer(end_chr1)
start_chr2 <- as.integer(start_chr2)

knownInsertion2 = GRanges(IRanges(c(start_chr1),c(end_chr1)),seqnames="chr1", chrB="chr2", startB=c(start_chr2), copied=TRUE)
knownInsertion2

sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
sim
metadata(sim)

Gives as an result:
> knownInsertion2
GRanges object with 4 ranges and 3 metadata columns:
      seqnames    ranges strand |        chrB    startB    copied
         <Rle> <IRanges>  <Rle> | <character> <integer> <logical>
  [1]     chr1  [21, 24]      * |        chr2        30      TRUE
  [2]     chr1  [16, 19]      * |        chr2        21      TRUE
  [3]     chr1  [ 6, 11]      * |        chr2        10      TRUE
  [4]     chr1  [ 9, 14]      * |        chr2        27      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
>
> sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                chr1
[2]    60 GGGGGGGGGAAAAAAGGGGGGGGGGGAAAACCCCCCAAAAAACCCTTTTCCCCCCCCCCC                                                                  chr2
> metadata(sim)
$insertions
  Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1    1 chr1     21   24 chr2     30   33    4   TRUE                 CCTT          TTCC
2    2 chr1     16   19 chr2     21   24    4   TRUE                 GGAA          AACC
3    3 chr1      6   11 chr2     10   15    6   TRUE                 GGAA          AAGG
4    4 chr1      9   14 chr2     27   32    6   TRUE                 CCAA          AACC

To make this customizable change:
sizes <- c(whatever, numbers, you, like)
output=NA (to a file path prefix to get a .csv file with the metadata)
genome = DNAStringSet() To replace the sample sequence
bpSeqSize = number (to change how many of the flanking region is shown in the metadata)


text:
Is there an R package that computes homoplasy excess ratios (HER)?
I want to calculate the homoplasy excess ratio for a phylogenetic tree, for example homoplasy resulting from recombination.
I am wondering if there are any R packages out there that will do this?
I have looked at phyloseq, Rphylopars, phylocurve, PhylogeneticEM, PhyloinformR, but these packages do not have functions that calculate HER. I searched "homoplasy excess ratio" on bioconductor, but it returned ca. 3000 packages.

Maynard's Smith's homoplasy is what you want and it is available in the recombination detection package (RDP) described here. Its really a Windows package rather than R.


text:
How to interpret duplication from MultiQC/FastQC?

FastQC of my sample files, aggregated into a single plot by MultiQC. Blue represents unique reads. Black represents duplicate reads. The x-axis is the number of reads.
I see I have more duplicate reads than unique reads, how do I interpret this plot? This is my attempt:

FastQC took the first 100,000 sequences in each sample (details).
  Generally, higher duplicated read percentage is welcome as it
  represents better sequencing depth and lower sequencing error
  contamination. This dsitribution in this plot looks reasonable and
  nothing unusual.

Can you provide us with some more info, including the typical number of reads in each library, the library type, source material, et cetera? It's hard to tell how to interpret your data without that info. If it's single-cell versus a sample from bulk tissue that makes a big difference in complexity. Also, something like RNA seq will have more duplicates that genome sequencing library types.

I would agree with your statement.  Typically the ratio of unique to duplicated reads is a function of input library complexity (number of input molecules) verse sequencing depth (the amount of molecules sequenced).
In the example you gave you have about a 2:1 ratio of duplicates to unique reads which suggests that these libraries have been sequenced to their fullest and additional sequencing would yield few new unique reads.

The higher number of duplicates could be in a high-complexity library sequenced very deep or in a low-complexity library sequenced with many fewer reads. Without more info from OP it is hard to interpret.

the x-axis provides the total number of reads.  either way it is the ratio that is important to tell you if the sequencing has saturated the library.


text:
Comparing between different set of clusters produced from the same network
I have a set of eight gene clusters which I obtained from clustering a network. However, I had to update the genes in the network at a later stage which lead to certain deletions (genes not present in the new gene list). I did the clustering again and this time there were a total of 9 clusters. Is there some way in which I could compare the earlier eight clusters with the new nine ones to effectively say that "clustering has happened along the same lines". It would have been easier to do an eight to eight comparison checking from the largest to the smallest cluster. However, since I have nine clusters now, I am kinda in a fix as how to compare.
I used the same clustering algorithm to cluster the network in both the cases

Answer from @llrs, converted from comment:
Perhaps I'm not understanding correctly what you are doing. But these stability tests are used when the number of clusters is unknown, so you test the resolution, or other metrics for a range of numbers of clusters. Then you choose the best one (either the one that is better by the score, or the one that increasing the number of clusters don't improve your metric much more)

If you have both cluster assignments, you can check for consistency of assignment by looking at pairwise counts using table in R:
cons.mat <- table(clusters.1, clusters.2)

The result of this will be a matrix showing the number of cells / samples that have the observed cluster assignments via each method. If the two cluster assignments are similar, then samples from one cluster via one method should be mostly assigned to one cluster in the other method. You can check this by comparing table cell values with their row/column sum:
## show consistency along the first dimension
apply(cons.mat, 1, function(x){max(x / sum(x))})
## show consistency along the second dimension
apply(cons.mat, 2, function(x){max(x / sum(x))})


text:
In calculating the retention index, why do we use the character state with the lowest frequency?
To calculate the retention index for a phylogenetic tree, we use the following formula:
$$\frac{\text{maximum number of steps on tree - number of steps on the tree}}
{\text{maximum number of steps on tree - minimum number of steps in the data}}$$
To calculate the maximum number of steps, we look at the number of observed character states for a character and take the lowest number. For each character, we take the lowest number of observed traits and then sum across all characters.
Take the following dataset (columns are characters, taxa are rows):
          123456
Outgroup  000000
Taxa A    100001
Taxa B    110001
Taxa C    111001
Taxa D    111101
Taxa E    111111
Taxa F    111110

The maximum number of steps for the first character is 1, since the character state 0 is attested once and character state 1 occurs six times.
What I don't understand is why we do not use 6 for the maximum number of changes for the first character.

Could you explain what is the matrix? Each column is a feature, but what does 0 or 1 represent, a step? Perhaps I would understand better if you could show the retention index of this dataset

The retention index is used to gain insight into how well the data conform to the most parsimonious tree. 0 and 1 represent the character state.

Agreed RI is used for a parsimony tree. The minimum of steps requires searching the total tree space of the data and is not trivial to obtain without an algorithm. The maximum number of steps would computationally derrived.

found this very useful. I was wondering you could help @Kohl, how do I go about working out max steps if I have a matrix with not only 0s and 1s, but other integers such as: eg taxonA 11000203044.
taxonB 11001102034 Thanks for your help, hope you see this after it being closed for so long

The maximum number of steps (or changes) is the number of taxa with state 1 or 0, whichever is smaller. Conceptually, this is the number of character changes if each taxon evolved its state independently of the other taxa.
If for a given character the character state 0 occurred twice and character state 1 occurred five times, 2 would be the maximum number of changes.

Thank you for this. Could you just edit your answer to explain why in the scenario you sketch above the five-time occurrence of character state 1 is not used in calculating the maximum number of steps? What I find confusing is that if character state 1 occurs five times, could that not represent five independent (i.e., homoplasious) character state changes?


text:
How to create a dendrogram of clusters and reconstruct phylogenetic relationship
I'm working with a Daphnia data set looking at the 16s gene pulled from the BOLD database. So far I did

a multiple sequence alignment and
attempted to cluster and plot a dendrogram:
myClustalWAlignment <- DNAStringSet(msa(MySeqString, cluster="nj", gapOpening=30.0,  gapExtension=20.0, maxiters=2, substitutionMatrix="default", type="default", verbose=TRUE))
Replace full sequence names with Accession numbers for easy identification
names(myClustalWAlignment) <- Daphnia_DF$Accession_number
Converts the alignment object into a DNAbin class. Some of the results are shown and may be useful for further analysis:
ClustalAlignment_dnaBin <- as.DNAbin(myClustalWAlignment)
Creating a distance matrix:
distmat <- dist.dna(ClustalAlignment_dnaBin, model = "TN93", as.matrix = TRUE, pairwise.deletion = TRUE)
clustering the sequences:
Daphnia_Clustered <- IdClusters(distmat,
                      method = "UPGMA",
                      cutoff = 0.01,
                      showPlot = TRUE,
                      type = "both",
                      myXStringSet = NULL,
                      model = MODELS,
                      collapse = 0,
                      processors = NULL,
                      verbose = TRUE)
checking to see the number of clusters:
length(unique(unlist(Daphnia_Clustered1)))

output: 80
The problem is that I want a dendrogram showing just the 80 clusters if that's possible. Right now the output of the above clustering is:

Is there a way to get a smaller number of clusters? Or even being able to plot just the 80 OTU clusters would be great. Is there an easy way of reconstructing the phylogenetic relationship between these clusters?

Welcome! Just curious - why don't you want to just make a maximum likelihood or Bayesian tree? Those are very standard approaches and you can just collapse nodes later on to make clusters.

Hey, I'm just extremely new r and bioinformatics I guess. I did attempt that and I did come out with a tree structure. However, how would I go on about collapsing the branches, are there any packages you particularly recommend ?

If you need to quickly visualize the tree I would just save the tree file, open it in FigTree, and collapse the nodes. https://www.biostars.org/p/18302/  If you add some info to your question like what the end goals are for this project or this phylogenetic analysis/visualization I think people can get a better answer tailored to the needs of the question.

I agree that it is easier to make a subtree in figtree and write the tree in R to the hard-drive.


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

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

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

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

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

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

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


text:
How to deal with duplicate genes having different expression values?
I have RNA-Seq data which is FPKM. In the dataframe df first column is gene_name and the other 100 columns are samples.
Usually if it is counts data I do like following:
df2 <- aggregate(. ~ gene_name, data = df, max)

I'm not sure what do with the FPKM data if there are duplicate genes with different FPKM value for the same sample.
Lets say it looks like this:
gene_name     sample1        sample2        sample3
5S_rRNA      0.3206147    0.3327312      0.377578
5S_rRNA      0.3342000    0.0000000      0.1305166

Any suggestions please.

How were the read mapped? Do you have access to the raw counts data?

I have downloaded TCGA FPKM data. Would like to use the data for co-expression analysis with cemitool

I assume you're familiar with the various issues surrounding FPKMs, so I'll not expound upon them.
As a general rule, you should be using gene IDs rather than gene names, since the former are unique while the latter are not. If you only have access to data quantified on gene names, then the appropriate way to merge RPKMs is with a weighted sum:
$FPKM_{gene} = \frac{FPKM_{copy1} * Length_{copy1} + FPKM_{copy2} * Length_{copy2}}{Length_{copy1} + Length_{copy2}}$
As an aside, rRNA expression levels will tend to be wrong, since they're normally excluded by either poly-A enrichment or the use of ribo-zero.

Sure I will use gene ids. but lets say I have counts data and before doing differential analysis do you think this is the right to way to filter duplicates; aggregate(. ~ gene_name, data = df, max)

If you have counts then sum them.

So, if the first column is Ensembl ids and all other columns are samples in a dataframe "df", Is this right way to sum the duplicate Ensembl ids? aggregate(. ~Ensembl_ids, data=df, FUN=sum)

There should never be duplicate Ensembl IDs, they're unique. You only have issues like this with gene names or UCSC IDs.

oh yes, so with gene names Is this the way to sum the counts of duplicates aggregate(. ~gene_names, data=df, FUN=sum)

I have to say I disagree with the weighted approach - why would you place more emphasis on any of the individdual loci (in this case the longer)?

The best way to deal with this is to use unique gene IDs, for example ensembl accession numbers. So use the ensemble gtf annotation when quantifying the read counts and not the gene symbols. Just to illustrate, when I look for "5S_rRNA" in ensembl's annotation, i see 18 different "genes" with that gene symbol. But which 2 you have is unclear now.
grep "5S_rRNA" ensembl_symbol.txt
"ENSG00000252830"       "5S_rRNA"
"ENSG00000276442"       "5S_rRNA"
"ENSG00000274408"       "5S_rRNA"
"ENSG00000274059"       "5S_rRNA"
"ENSG00000276861"       "5S_rRNA"
"ENSG00000274759"       "5S_rRNA"
"ENSG00000280646"       "5S_rRNA"
"ENSG00000277411"       "5S_rRNA"
"ENSG00000201285"       "5S_rRNA"
"ENSG00000212595"       "5S_rRNA"
"ENSG00000277418"       "5S_rRNA"
"ENSG00000277049"       "5S_rRNA"
"ENSG00000274097"       "5S_rRNA"
"ENSG00000277488"       "5S_rRNA"
"ENSG00000274663"       "5S_rRNA"
"ENSG00000283433"       "5S_rRNA"
"ENSG00000275305"       "5S_rRNA"
"ENSG00000278457"       "5S_rRNA"

Yes you are right. But can I take median of those duplicates? something like aggregate(. ~ gene_name, data = df, median)

Expanding this answer, it seems that it was mapped to ensembl ids. So I would recommend to keep them as they are, and only when looking for biological meaning to change to gene symbols (or when talking to others)

I would strongly advice not to use FPKM at all, it is pretty clear now that they are not the right values to use for any analysis. If you still want to use it, keep your two genes apart like they were.


text:
Bioconductor, how to retrieve the genes from a GEO gene expression dataset?
I'm quite new to Bioconductor and GEO.
Here's my task: I've to develop a simple R script able to read a gene expression dataset on GEO and extract the genes related to it. But... I don't know the command to extract the gene list.
The GEO dataset I want to analyze is GSE59867.
Here's my R code:
options(stringsAsFactors = FALSE)

library(oligo)
library(GEOquery)
source("https://bioconductor.org/biocLite.R")
biocLite("affyio")
library("affyio")

getGEOSuppFiles("GSE59867")
list.files("GSE59867")

untar("GSE59867/GSE59867_RAW.tar", exdir = "GSE59867/CEL")
list.files("GSE59867/CEL")

celfiles <- list.files("GSE59867", full = TRUE)
rawData <- read.celfiles(celfiles)

Now I have the rawData containing the gene expressions.
How can I extract the list of genes associated to rawData at this point?
Thanks!

Have a look at the pd.hugene.1.0.st.v1 package. That contains the annotations for that microarray.


text:
How to transform a DNAStringSet from the Bioconductor package Biostrings to a data frame?
I am working on Mac OS X. I am using R version 3.5.1.
I have important a FASTA file into R using Biostrings::readDNAStringSet. This creates a DNAStringSet object in R.
I would like to transform it into a data frame so that I can then plot the widths with ggplot2.
as.data.frame on my DNAStringSet object returns a data frame with only one column. Is there a function that returns a data frame with columns width, seq and names present in the DNAStringSet object?

According to the documentation (?Biostrings::DNAStringSet):

width(x): A vector of non-negative integers containing the number of
            letters for each element in x.  Note that width(x) is
            also defined for a character vector with no NAs and is
            equivalent to nchar(x, type="bytes").
names(x): NULL or a character vector of the same length as x
            containing a short user-provided description or comment for
            each element in x.  These are the only data in an
            XStringSet object that can safely be changed by the user. All
            the other data are immutable!  As a general recommendation,
            the user should never try to modify an object by accessing
            its slots directly.
as.character(x, use.names=TRUE): Converts x to a character vector
            of the same length as x.  The use.names argument controls
            whether or not names(x) should be propagated to the names
            of the returned vector.

This should be enough for you to create a data.frame:
dss2df <- function(dss) data.frame(width=width(dss), seq=as.character(dss), names=names(dss))

Also, if you have a huge DNAStringSet, perhaps you do not want to convert all the sequences to characters; storing sequences as characters supposedly is less efficient. For plotting it might be sufficient to have just widths.
qplot(width(dss), geom='histogram'))


text:
Error in quantile.default(x, p = p) in EgdeR calcNormFactors
I am trying to run TMM normalization using rpy2 and when I run calcNormFactors() function:
dge_list = edgeRLib.DGEList(counts=rawcounts)
dge_list = edgeRLib.calcNormFactors(dge_list, method="TMM")

I am getting the error:

RRuntimeError: Error in quantile.default(x, p = p) :
    missing values and NaN's not allowed if 'na.rm' is FALSE

So, I thought that maybe my rawcounts 2D numpy array contains nan values but it does not. I checked it with both pd.isnull() and np.isnan() functions. So, I suspect now that the 2D array itself can not be normalized for some reason and I do not know why. The error is most probably not related to rpy2 set up and R installation because it was working before for similar numpy arrays but failed this time. Any suggestions would be greatly appreciated.

Ok, it turned out that some rows and/or columns contained all values as zeros and that was the reason for the error. I wrote a function to fix that:
def fix_rawcounts(rawcounts, cell_ids, gene_symbols):
    '''The function checks whether there are
    any rows or columns with all zeros
    and removes them'''
    print('Checking whether rows or columns contain all zeros')

    # Checking columns
    count_arr = np.apply_along_axis(np.count_nonzero, axis=0, arr=rawcounts )
    rawcounts = rawcounts[:,count_arr != 0]

    # Checking rows
    count_arr = np.apply_along_axis(np.count_nonzero, axis=1, arr=rawcounts )
    rawcounts = rawcounts[count_arr != 0,:]

    print('finished checking zero rows and columns')
    return rawcounts


text:
Bash scripting FastQC for multiple fastq files in multiple directories
I am completely new to bioinformatics so I'm looking to learn how to do this.
I have multiple directories with fastq files: E.g; 10 Directories with each time series, each with Treatment and control directories, each with rep1 rep2 rep3.
For example: T9/Infected/Rep1/*.fastq.gz.
I'm looking to create a loop to run fastQC on each fastq file instead of having to submit a separate job for each directory.
Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.

multiqc kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the multiqc report and the actual fastqc html report to get a better idea of any error modes.
Going off of @Kubator's answer, I noticed that there was no command to run fastqc.
Here's a simple one-liner to run fastqc in parallel on all of your fastq files. The -j 25 uses 25 threads. Change 25 to however many threads you want/have for max speed.
#Run fastqc on everything in parallel.
> find ../reads/ -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' | parallel -j 25 --verbose
# copies all the fastqc files to directory ./
> find ../reads/ -name '*fastqc.*' | xargs -I '{}' mv '{}' ./

These files might be output by multiqc, anyway!
edit - now the command is safe for whitespace in file names. Please, no newlines in filenames though!

+1, but note that this will break in the unlikely case where any of the fastq file names contain whitespace. A safer version is find . -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' but that still fails in the (even more unlikely) case where a file name contains a newline. This should work for anything (but requires a version of find with -printf, like GNU find): find . -name '*.fastq.gz' -printf '"%p"\n' | parallel -j 25 --verbose.

whitespaces in file names. Oh, the horror! Thanks for your comments!

heh, I know. Anyone who has spaces in fastq filenames deserves what they get :) sorry, I'm just very used to posting on [Unix.se] where this sort of detail is more relevant.

This is the typical kind of job for snakemake.
Assuming your have one file per replicate named for instance T9/Infected/Rep1/Rep1.fastq.gz, you can prepare a file that you call Snakefile with the following content:
timepoints = list(range(10))
conditions = ["control", "infected"]
replicates = [1, 2, 3]

rule all:
    input:
        expand(
            "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html",
            time=timepoints,
            cond=conditions,
            rep=replicates)

rule do_fastqc:
    input:
        fastq = "T{time}/{cond}/Rep{rep}/Rep{rep}.fastq.gz"
    output:
        html = "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html"
    shell:
        """
        fastqc {input.fastq}
        """

Put this file in the directory that contains the T* directories and run snakemake from there.
The top all rule explains which files you want. The do_fastqc rule explains how to make one fastqc report from one fastq.gz file.
With a bit more work, this can be used to submit jobs to a computing cluster. Snakemake has some tools for this.
If you don't know the exact names of the fastq files but they all follow the same pattern, you will need to use the glob python module and do a little bit of programming to determine the possible values for rep, cond and time. The "snakefile" can contain any python code you want.
If there are no regular pattern in the file names, fix this issue first ;)

Example dir structure:
$ find FastQC/
FastQC/
FastQC/T9
FastQC/T9/Infected
FastQC/T9/Infected/Rep1
FastQC/T9/Infected/Rep1/test11.fastq.gz
FastQC/T9/Infected/Rep1/test1.fastq.gz
FastQC/T9/Infected/Rep2
FastQC/T9/Infected/Rep2/test2.fastq.gz
FastQC/T9/Infected/Rep3
FastQC/T9/Infected/Rep3/test3.fastq.gz

If understood well You need to run some job on every *.fastq.gz file. Then You can do something like this (my example job is gzip test, replace with your job):
Rookie:
$ find FastQC/ -type f -name "*.fastq.gz" | xargs gzip -tv
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

Solid:
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} gzip -tv {}
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

find finds files with name *.fastq.gz and outputs it with zero byte delimited (to support weird characters like space etc. in filenames)
xargs represents output as {} and passes it to gzip -tv

If You want to copy files inside one heap folder:
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} cp -pv {} FastQC_heap/
`FastQC/T9/Infected/Rep1/test11.fastq.gz' -> `FastQC_heap/test11.fastq.gz'
`FastQC/T9/Infected/Rep1/test1.fastq.gz' -> `FastQC_heap/test1.fastq.gz'
`FastQC/T9/Infected/Rep2/test2.fastq.gz' -> `FastQC_heap/test2.fastq.gz'
`FastQC/T9/Infected/Rep3/test3.fastq.gz' -> `FastQC_heap/test3.fastq.gz'

You can use following command from the parent folder to run fastqc on fastq.gz files in subdirectories
find -name '*.gz' | xargs fastqc

what is the reason for downvote? Can you mention the error in the command

I think you could also explain what it actually does (but I did not downvoted you, don't know)

dear downvoters, this the best method that can be used to run fastqc on the files located in different locations...

Also, what's the advantage of your approach compared to https://bioinformatics.stackexchange.com/a/5295/57?

@KamilSJaron are you asking me the advantage of this approach?. This is an alternate way of doing it. This is quite faster compared to other methods and command is also easy. The requirement of the person here he mentioned as Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.


text:
Simulating 3' end tag-based scRNA-seq reads
Are there any tools that will simulate 3' end tag-based single-cell RNA-seq reads? For example, Drop-seq, 10X Chromium, CEL-seq?
There are tools that simulate scRNA-seq gene count data (e.g. Splatter), but I can't find anything that will simulate reads or even counts at the transcript level.

Looks like a couple options have become available:
Minnow: article code
scReadSim: preprint code


text:
ScaleData() from Seurat causes crash on RStudio Cloud
Using the sample data in the 2,700 PBMC clustering tutorial, the session crashes at the ScaleData() step. The progress bar remains at 0% and hangs, followed by an R crash. I've tried reducing the size for number of genes to scale at in a single computation with the argument block.size with no change. The Seurat S4 object is only ~70mb so I can't imagine I'm exceeding the RStudio Cloud 1gb RAM limitations.
If you'd like to check out the cloud instance I can add you to the workspace temporarily. There is a limit on members allowed in a workspace.
I've been trying to run this tutorial on another remote cluster in RStudio inside of a container (Singularity here vs. Docker in RStudio Cloud) and am hitting the same snag.
I get receive no errors from RStudio Cloud, but in the other evirnment I get this error:
ERROR session hadabend; LOGGED FROM: rstudio::core::Error {anonymous}::rInit(const rstudio::r::session::RInitInfo&) /home/ubuntu/rstudio/src/cpp/session/SessionMain.cpp:563.
On this issue thread someone mentions it they receive this error message when hitting their RAM limits.
The source for ScaleData() can be found in the Github repo.

Are you following the tutorial exactly? I don't see ScaleData in the tutorial.

RegressOut() is a deprecated function, regression now happens in the data scaling step. vars.to.regress is now passed as an argument to ScaleData(). Somehow you've got an old version of the tutorial, see this tutorial compiled 22 March, 2018.

It seems RStudio server doesn't know it's limitations. Currently, Projects are limited to 1GB of RAM. Using the R package ulimit which mirrors the function of ulimit and allows implementation of memory limits in linux I was able to set my own limitations within the R session. After setting the limit to 1GB I see the ulimit limitation being tripped when invoking ScaleData(). Looks like a case of exceeding available system resources.


text:
hgsql not found when running a script despite being accessible on command-line
I am working on macOS. I am following the steps described here to lift over an annotation from the reference genome of C. elegans to another assembly of C. elegans.
I have run this code:
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

However, the do.log file contains:
zsh:1: command not found: hgsql
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Sun Nov 11 15:38:41 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/job.csh
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
# ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
cd /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat
rm -rf tParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 58018742 0 /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.chrom.sizes 2000 -lstDir=tParts
rm -rf qParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 10000000 0 ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes 1000 -lstDir=qParts
Can't open seq size file ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes: No such file or directory
Command failed:
ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh

I think the biggest problem is that the hgsql command is not found.
If I run hgsql on the command line, I find it even though I get the following error:
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
  Referenced from: /usr/local/bin/hgsql
  Reason: image not found
[1]    34059 abort      hgsql

Location of the hgsql installation:
which hgsql
/Users/username/Documents/projects/main/output/software/liftover/data/bin/hgsql

Can anyone help me with using hgsql, please?

It is possible that the shell being used in the process is different from your standard command line shell. I am not familiar with the tool, and it's hard to say more without knowing your OSX version, but it looks like it is using zsh. I think that zsh only fairly recently became the standard OSX shell. I don't know the tool or the shell very well, but it looks like you might be able to explicitly set up hgsql in your PATH in your .zshrc.
It's also possible that the liftover script was not built to run via zsh. I notice that it is trying to run csh scripts as part of the workflow. Then again the localhost may insulate it.
So I think that my first step would be to sort out which shell I am using, whether it works with the liftover pipeline, and how it's running on the localhost. I notice that the ssh call running the csh script on the localhost is the step that actually fails. What shell is executing the csh script? Does that shell know about hgsql? It could be that the way that the error bubbles up it is hiding which shell it is actually using.
I'd recommend trying that ssh call in isolation from the rest of the pipeline, and diagnosing the error that way.
Hope that helps.


text:
Continuing analysis based on a log-likelihood ratio? (Comparing methylation between two PoolSeq phenotypes)
I'm trying to get some pointers on where to start on the analysis of my dataset. First, some basic info:
Background:
We performed Nanopore sequencing of two poolseq samples called UU and DD referring to damage vs undamaged. We were able to obtain approximately 3x coverage (this is currently more of a proof of concept than a hardcore association job, so the coverage is unfortunate but worth pursuing). Between these samples, we know that one is going to have different methylation patterns in response to a phenotypic trait response, but we don't know what regions along the genome are going to be differently methylated. The goal is to compare what's methylated (as in higher or lower levels of methylation) between the two samples, knowing there ought to be a lot that cancels out and come up with ideally a short list of significantly different regions we can associate to the plastic response.
Output of initial methylation calling data:
Using the raw signal level data against the reference genome and sorted, mapped BAM files the Nanopolish software was used to call methylation. This produces a dataset with the following headers:
1. Chromosome, 2. +/- Strand 3. Start vs. 4. End of the CG dinucleotide (the difference between these is 1 or two bases) 5. The read ID the dinucleotide came from 6. the log likelihood ratio calculated from the model embedded in the methylation calling process using 7. log likelihood that it's methylated and 8. log likelihood that it's unmethylated. 8. number of calling strands (always appears to be 1) 9. number of motifs (don't know about this) and finally 10. the sequence surrounding the dinucleotide, which is always short but variable in length.
I can follow the basic interpretation of the log-likelihood... if it's a positive number, there is statistical support that it's methylated. This would probably be about as much as we need, short of a p-value, if we were looking for a specific region, but we are trying to analyze the differences in de/methylation across the entire methylome. How would you suggest using the statistic values provided in our dataset to compare and visualize the significant different modifications of the dinucleotides between strands/samples?

I assume you don't want to spend 1-2 years developing your own method based on nanopore data, so the path forward is to convert the output into a format used by standard DMR calling programs. There's no standard format for that, annoyingly, but the general things you'll need are:

Coordinates of a C/CpG with coverage
Number of reads supporting methylation at that site.
Number of reads not supporting methylation at that site.

Typically you would do some filtering before/after making such a table so you only have decent quality calls (that would translate into some reasonable log likelihood for your data) and possibly only sites with decent (5x, 10x, 30x, etc.) coverage. You'd then feed that for all of the replicates in all of your groups into programs like bsseq, DSS or Metilene (to not so randomly select 3 examples).
If there aren't DMR finding programs already developed specifically for nanopore data I expect someone is developing it. The technology has a higher error rate, but given that bisulfite conversion as required in Illumina data also leads to a higher error rate that's not necessarily a huge issue if the methylation calling itself is reasonably accurate and the depth is high enough. Note that you probably need 10X coverage as a bare minimum unless your effect size is very very large.


text:
In Uniprot how do I know which variants are mutations and which are polymorphic?
In Uniprot how do I know which variants are mutations and which are polymorphic?
I am looking at P38398 (BRCA1_HUMAN) specifically.
Many thanks.

What do you mean by "mutation"? Are you just talking about how frequent a variant is in a population?

This answer is very late but may still help others:
Have a look at the UniProt index of human polymorphisms and disease mutations at https://www.uniprot.org/docs/humsavar
Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.

Thanks for answering! No answer is ever late here, the objetive of this site to to amass a collection of useful Q&As, so the age of a question is irrelevant as long as it could still be asked by someone else. If, as your answer seems to suggest, you are affiliated with UniProt, it would be helpful if you could mention that since it would add more weight to your answers.


text:
E.coli Sequencing & Analysis
I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.
The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.
I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.
For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.
Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!
Thanks,
Alex.

You can first do assembly of contigs with spades and then use a reference to get longer sequences, for example you can try this: https://github.com/fenderglass/Ragout

A de novo assembly with spades should work as well as a mapping with bowtie2.
While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.
When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.
I'd like to give you the advice to use some kind of distribution system like Anaconda which helps you to keep your system clean and make it easier to reproduce your steps.
Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :
De novo assembly
1. Quality Check ( you already did this )
2. Run Spades
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log

This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like.
You now have a folder called k12 which holds the assembly information.
To get some information about the contigs you can use the emboss toolkit.
infoseq k12/contigs.fasta

fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 -              NODE_859_length_78_cov_1.000000 -              N    78     43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 -              NODE_860_length_78_cov_1.000000 -              N    78     44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 -              NODE_861_length_78_cov_1.000000 -              N    78     43.59

This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.
3. Remove low quality contigs
There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.
4. Check for possible contamination
So to check for possible contamination you can use kraken2. But it needs some disk space. I am using here the minikraken database
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log

The ecoli_kraken.report shows an overview of how the contigs could be classified to availible species from the database.
Based on this results you can also remove some contigs.
5. Run Prokka
Now after you removed unwanted contigs, you can run prokka to get the annotation.
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log

And you will find the information in the prokka folder
Mapping
The mapping process is a little bit easier. You can use bowtie2 and samtools
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam

158216 reads; of these:
  158216 (100.00%) were paired; of these:
    35696 (22.56%) aligned concordantly 0 times
    120047 (75.88%) aligned concordantly exactly 1 time
    2473 (1.56%) aligned concordantly >1 times
    ----
    35696 pairs aligned concordantly 0 times; of these:
      4458 (12.49%) aligned discordantly 1 time
    ----
    31238 pairs aligned 0 times concordantly or discordantly; of these:
      62476 mates make up the pairs; of these:
        56512 (90.45%) aligned 0 times
        5365 (8.59%) aligned exactly 1 time
        599 (0.96%) aligned >1 times
82.14% overall alignment rate

You now also get some information about the similarity to your reference sequence.
Now you need to work a little bit more on the BAM file. You need to sort and index this.
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam

Now you can create a consensus sequence using bcftools and seqtk
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta

This gives you now a sequence where unknown bases are marked as N.
As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file bio.yaml
name: bio
channels:
  - bioconda
  - defaults
dependencies:
  - aragorn=1.2.38=h470a237_2
  - barrnap=0.9=2
  - bcftools=1.9=h4da6232_0
  - bedtools=2.27.1=he941832_2
  - blast=2.6.0=boost1.60_0
  - bowtie2=2.3.4.3=py36h2d50403_0
  - emboss=6.6.0=h6debe1e_0
  - hmmer=3.2.1=hfc679d8_0
  - infernal=1.1.2=h470a237_1
  - kraken2=2.0.7_beta=pl526h2d50403_0
  - libdeflate=1.0=h470a237_0
  - libidn=7.45.0=2
  - minced=0.3.2=0
  - parallel=20160622=1
  - perl-app-cpanminus=1.7044=pl526_1
  - perl-bioperl=1.6.924=4
  - perl-carp=1.38=pl526_1
  - perl-constant=1.33=pl526_1
  - perl-data-dumper=2.161=pl526_2
  - perl-encode=2.88=pl526_1
  - perl-exporter=5.72=pl526_1
  - perl-extutils-makemaker=7.34=pl526_2
  - perl-file-path=2.15=pl526_0
  - perl-file-temp=0.2304=pl526_2
  - perl-parent=0.236=pl526_1
  - perl-threaded=5.22.0=13
  - perl-xml-namespacesupport=1.12=pl526_0
  - perl-xml-parser=2.44=pl526h3a4f0e9_6
  - perl-xml-sax=1.00=pl526_0
  - perl-xml-sax-base=1.09=pl526_0
  - perl-xml-sax-expat=0.51=pl526_2
  - perl-xml-simple=2.25=pl526_0
  - perl-yaml=1.27=pl526_0
  - prodigal=2.6.3=1
  - prokka=1.12=3
  - samtools=1.9=h8ee4bcc_1
  - spades=3.12.0=1
  - tbl2asn=25.3=0
  - boost=1.60.0=py36_0
  - bzip2=1.0.6=h14c3975_5
  - ca-certificates=2018.03.07=0
  - certifi=2018.10.15=py36_0
  - curl=7.62.0=hbc83047_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - giflib=5.1.4=h14c3975_1
  - icu=58.2=h9c2bf20_1
  - jpeg=9b=h024ee3a_2
  - libcurl=7.62.0=h20c2e04_0
  - libedit=3.1.20170329=h6b74fdf_2
  - libffi=3.2.1=hd88cf55_4
  - libgcc=7.2.0=h69d50b8_2
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgd=2.2.5=hceca4fd_3
  - libpng=1.6.35=hbc83047_0
  - libssh2=1.8.0=h1ba5d50_4
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libtiff=4.0.9=he85c1e1_2
  - libuuid=1.0.3=h1bed415_2
  - libwebp=1.0.0=h222930b_1
  - libxml2=2.9.8=h26e45fe_1
  - ncurses=6.1=hf484d3e_0
  - openjdk=8.0.152=h46b5887_1
  - openssl=1.1.1=h7b6447c_0
  - perl=5.26.2=h14c3975_0
  - pip=18.1=py36_0
  - python=3.6.7=h0371630_0
  - readline=7.0=h7b6447c_5
  - setuptools=40.5.0=py36_0
  - sqlite=3.25.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - wheel=0.32.2=py36_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=ha838bed_2

Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.
conda env create -f bio.yaml

I hope this gives you some help how to process your data. If you need some information just let me know ;)


text:
Double-counting coverage of overlapped read pairs
EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.
I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent).
I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:

Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?
samtools (v1.7):
samtools depth -a -b exome_targets_GRCh37.bed sample1.bam > sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt

There is already code in htslib for resolving overlapping paired-end reads, but it only seems to be used in samtools mpileup and not in samtools depth. I guess you could modify bam2depth.c and add the bam_mplp_init_overlaps() function call. I used it in my own program for strand-specific coverage calculation, see here

That is interesting indeed and I will definitely look into the bam2depth.c modifications in order to drop overlap from the function call. Was it difficult to implement this into your program? I see your program and the overlapping ends counting feature appears to be what I am looking for. I will need to locally install samtools as I work on a cluster. Thanks for the recommendation!

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use samtools mpileup instead:
$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

I haven't shown the entirety of the output and you'd want to use the -aa flag to ensure that all bases are returned.
Note that samtools mpileup is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to -Q 1 (the default is -Q 13, which you'd want to change).
As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use bamPEFragmentSize from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.

Does this method take bases that are not covered by a read into account?

Yes, that's what the -aa flag I mentioned does.

I tested this on data and it seems to handle the issue nicely! Thanks Devon.

One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique sequence.
Programs like SeqPrep, PEAR (Paired-End reAd mergeR), and fastq-join can do this fairly quickly. However, this isn't a perfect solution as the reads will not be joined unless they meet some minimum overlap criteria.

added an edit, I do not want to remap reads or make any modifications to existing reads.

picard CollectHsMetrics has the parameter PER_BASE_COVERAGE which take this into account.
You can run it like this:
 $ java -jar picard.jar CollectHsMetrics \
      I=input.bam \
      O=hs_metrics.txt \
      R=reference_sequence.fasta \
      BAIT_INTERVALS=bait.interval_list \
      TARGET_INTERVALS=target.interval_list \
      PER_BASE_COVERAGE=base_coverage.csv

BAIT_INTERVALS and TARGET_INTERVALS could be the same file. You will need to convert your bed file to it using picard BedToIntervalList.
Another useful tool is mosdepth.
$ mosdepth --by regions.bed SampleName input.bam

fin swimmer

Forgive me because I’be only looked into the tool, never used it. Do I have to generate some kind of sample.dict file for every sample I want to do this on or is it straightforward. It’s no big deal doing this on a reference but it would have too much overhead running it on every sample. Could you provide the steps involved in running this?

A sample.dict isn't needed. I edit my answer as @terdon asks for.

Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream tools. For example, the merging tool of USEARCH increases base quality scores of bases that match between overlapping reads. Depending on the context under which coverage information will be used in the future, it may be a good idea to "double-count" overlapping bases.

The reads aren't from different clusters in the traditional sense, if that were the case they wouldn't be assigned as mates. The traditional method to exploit overlapping reads is to increase/decrease phred scores on one of the mates according to whether the base calls in overlapping regions are the same or different, respectively. That's performed internally in things like samtools mpileup.

I agree. The reads come from the same cluster.  However, you could match up mates from different clusters. The Illumina sequencer uses a fragment index to match up mates from the two mate sequencing runs. It also appears from your comment that samtools incorporates the extra information gained from having two reads cover the same base. Just so we're all on the same page, here is a video describing Illumina paired-end sequencing: https://youtu.be/fCd6B5HRaZ8


text:
Problem: "Pair-end" reads scRNA seq data (Drop-seq)
In case of Drop-seq, we have paired end data.
Read 1: Cell code + UMI (unique molecule identifier)
Read 2: The transcript information
But I have a problem/doubt with the sample I am working on.
The sample I am using is the following (Check the "Reads" tab).
As you know the Drop-seq is "paired-end", we are expected to see two reads for a spot. Although this sample say paired-end, it has only one read per spot.
For example I can share a link of a different scRNA-seq data where you can properly see two reads for a spot
Example sample, you need to check the "Reads" tab.
Where I am going wrong?
I asked one of the main authors of the paper. The following is the reply I got :

I recommend that you download the aligned BAM files that are hosted in
  the same GEO record. Read 1 is already processed into the cell and UMI
  barcodes and held as custom tags (XC and XM) in the BAM files. The
  cells are already barcode-corrected, so if you use those files, your
  cell barcodes will line up with mine; if you start from FASTQs, they
  will not. For most aligners, you can just use the BAM file as input to
  realign. (It has all reads, even those that did not align.)

But I could not find any "XM" or "XC" keywords in the bam file :(
To understand his reply you have to be familiar the processing steps of the Drop-seq.
Looks like they have submitted some kind of processed data. I could not figure out how much the data is processed. I am trying to use the data starting at different processing steps but I could not figure out how much the data is processed.

Where did you even find their BAM files in GEO? I only see the output from CellRanger (UMI counts and such).

My guess is the authors uploaded BAM files to GEO which have been converted back to fastq in SRA, losing all the extra annotation.

I received the following message from the author of another Drop-seq paper:

After much conversation with the folks at the SRA, it appears the SRA
  compression algorithm automatically deletes certain subsets of bam
  tags and does so in cryptic way (in the sense of not alerting the
  uploader). Obviously this is going to be an issue for many single cell
  datasets encoded this way, and I hope the SRA will eventually change
  their system to accommodate.

If you look at the "Download" tab of the same page, there now exists a link to download the originally uploaded BAM file.
Alternately, you can use the SRA AWS tutorial to download the files on that platform.


text:
Can I get the graph generated by cellranger
I’ve run the cellranger analysis pipeline on single cell RNASeq datasets. I can import the matrix and graph-based clusters into R. Doing this I can optimise the dimension reduction and plot cells with coloured by clusters generated by cellranger. I wish to optimise the graph-based clustering as well.
Is it possible to either:

obtain the graph generated by cellranger from the files generated. (I can only find the final clusters in the analysis output)
reproduce generating a similar graph (used for clustering) from the output matrix (preferably with R or Python.

With either approach, I am looking for a list of edges for nearest neighbors in the gene-barcode UMI matrix to use as input for a clustering algorithm.

Could you provide some data as an example or pick the data from the vignette to answer you would work for you? What have you tried to get the graph or reproduce it?

I’ve tried several ways to reproduce this graph in R but I’m unsure that I understand well enough what cellranger is doing (based on Python code in the GitHub repo) to generate inputs for the Louvian algorithm. I’m now unsure that this is necessary since any cellranger run generates various output files for the analysis in a specific directory structure. I think this question requires someone with familiar with this software already and has run it themselves. I cannot provide the data that I’m working with currently and the files are too large even if I could.

Hi @TomKelly. I don't know if I understood exactly the question. However, you may find useful to use scanpy for importing the data, run the Louvain clustering, and plot the data with the clusters. That is a Python module.

@TomKelly Could you provide your attempts and why they fail? This would make it easier for me to avoid this pitfalls (if I try to answer your question)

I don't understand what plot you want. You know where the T-SNE coordinates are, right?

No, I can perform tSNE in R. I want the graph that the graph-based clusters are derived from. I want to annotate clusters into it. Can these be found in the output data or does it have to be computed from the matrix?

Cell Ranger
You can't download the tSNE coordinates for cells directly from the Analysis tab of the fancy, polished .html document that Cell Ranger produces. If you have access to the machine on which the pipeline was run, you can grab the intermediate files. Everything is well documented by 10X. They explain defaults etc. on that page. The tSNE coordinates are available at:
analysis/tsne/2_components/projection.csv

and the contents look like:
$ head -5 analysis/tsne/2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300

Seurat
And the second point: Cell Ranger is replicating the Seurat workflow. Seurat uses the Louvain method for community detection by default. Check out the intro tutorial and see if you can follow along. This tutorial works with 10X data, in fact the first function call in the tutorial reads in the 10X matrices. If you are already plotting in R you'll speed through this tutorial. Seurat is great to work with, you'll have tSNEs and clusters in no time. And there's so much more you can do from there!

Yes I’m familiar with the Seurat workflow. I don’t understand the Python source code well enough to know which parameters were used by cellranger.

If you're unwilling to use a different software package, why mention that as a possibility in your second point? The matrices are very lightly analyzed data. The unfiltered matrices are counts for all cell barcodes. The filtered matrices are count for cell arcodes determined to come from real cells. You still have to normalize and perform quality control filtering on this data efore it's ready for its first round of dimensional reduction and clustering.

If you want information about their algorithms, see here. There is also info on that page about the matrices. This should help with your understanding of how Cell Ranger is massaging the data.

I’m considering whether it is necessary to compute the graph again of it has already been done by cellranger. If it has not been saved by this program, it will be needed to compute it as you’ve suggested. In that case, more detailed information will be helpful. I have (and will continue) to do it this way but it would be reassuring to know I’m on the right track and not reinventing a graph structure that I already have saved to disk somewhere.

Ah, I was thinking you just just had the final product. If you have access to the machine where the pipeline was run you can get the tSNE coordinates.

Yes I ran the cellranger pipeline myself. However, I’m not looking for the tSNE coordinates. I’m looking for a list of edges used for the graph-based clustering (I can already find the clusters).

I've managed to get want I needed with the Seurat::FindClusters function from the object@snn output. I've detailed my process below.

Files in a cellranger run
Cellranger removes the intermediary files after the clusters have been computed. This can be observed in the following path from the cellranger project directory. This contains empty directories.
SC_RNA_AGGREGATOR_CS/SC_RNA_ANALYZER/RUN_GRAPH_CLUSTERING

However, the clusters generated by cellranger can be found in the file:
outs/analysis/clustering/graphclust/clusters.csv

Clusters can also be "exported" from outs/cloupe.cloupe file into a CSV using the Cloupe Browser. These CSVs can be imported and used for plotting in R or Python.
The graph structure used to compute these (with the Louvian algorithm) cannot be found in the files generated by cellranger.
Reproducing graph-based clustering
The Suerat R package performs graph-based clustering. The shared nearest neighbour (SNN) modularity optimisation uses the same parameters as the Python scripts called by cellranger (k=30). This generates an SNN matrix which can be used to generate a graph structure from an adjacency matrix. For a Seurat object test the SNN sparse matrix can be extract with:
FindClusters(test, dims.use = 1:20, algorithm = 1,
             resolution = 0.6, k.param = 30, save.SNN = TRUE,
             reduction.use = "pca", metric = "euclidean")
SNN_Matrix <- as.matrix(ceiling(as.matrix(test@snn)
graph_structure <- igraph::graph_from_adjacency(SNN_matrix)

FindClusters from the Seurat package requires a Seurat object that has already has a dimension reduction run on it such as RunPCA. The algorithm = 1 is the Louvian algorithm as performed by cellranger. The algorithm, reduction.use and metric can be changed as mentioned in the Suerat package documentation. Clustering algorithms not yet supported by Seurat can be run on the graph_structure resulting from extracting the SNN matrix. It is not necessary to compute the edges separately, the Seurat functions will compute and return these.
For example, the Leiden algorithm can be computed as follows:
devtools:install_github("TomKellyGenetics/leiden")
library(leiden)
clusters <- leiden(SNN_Matrix)


text:
What is "block-compressed" file in samtools?
SAMtools returned an error message for my gzipped genome FASTA:

Indexed block-compressed FASTA file cannot be handled

The source code for the message is here.
// check if the it is a valid block-compressed file
if (IOUtil.isBlockCompressed(path, true)) {
    throw new SAMException("Indexed block-compressed FASTA file cannot be handled: " + path);
}

What's a "block-compressed" file?

A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block so that the whole file does not need to be decompressed and therefore can be accessed randomly. This isn't true of files compressed as a single block using gzip.
As far as I was aware samtools only works with uncompressed FASTA files.


text:
Using Seurat to compare mutant vs.wt
I am interested in using Seurat to compare wild type vs Mutant. I don't know how to use the package. How can I test whether mutant mice, that have deleted gene, cluster together?

This question is too vague and open-ended for anyone to give you specific help, right now. Can you provide a bit more detail? Have you run Seurat on your WT samples and have well-defined clusters that you want to compare against in your mutant samples?
If you can show some figures or briefly describe what you've done so far, or specifically what you're trying to achieve, it will be much easier for someone to provide guidance.

Single-cell analysis to compare samples is a long a difficult process. There is very good documentation for 10x Genomics cellranger, the DropSeq Pipeline and the Seurat R package. These tools all have GitHub repositories and the authors are very responsive if you encounter issues.
Depending on the technology used to generate the data, you'll need to use either cellranger or DropSeq to process the FASTQ files. These are designed to account for the the different experimental design. Still, it is important to ensure that your samples have been processed in the same way to reduce batch effects.
This question is to open-ended to cover all of the details needed to perform the analysis but here is some of the main things to consider and the general process to perform such as analysis. You will need competence using the command-line (shell) and programming in R (or Python) to perform single-cell analysis. I recommend performing this analysis on a remote server as some steps are memory intensive.
The general overview to compare samples is:

Demultiplex the reads for different indicies. You need to do this even if you've only put one sample per lane (to remove the index from the reads). cellranger mkfastq or Illumina's bcl2fastq will do this.
Obtain a reference genome (FASTA) and gene annotation (GTF) for the species you are working with. You can prepare a reference transcriptome with cellranger mkgtf and cellranger mkref.
Run cellranger count or the DropSeq pipeline on each sample separately. These will both perform STAR (splice-aware) alignment of paired-end RNA-Seq reads and count UMIs for each cell-barcode and gene. Cells will be identified and filtered by an automatic UMI threshold (this can be changed by forcing with expected number of cells). This will return a gene-barcode matrix.
You need to aggregate results of cellranger runs for different samples with cellranger aggr. This will perform downsampling by default to normalise the number of reads between samples. It will also perform graph-based (Louvian) clustering on the combined samples and return HTML and Cloupe summary data to explore, in addition to a combined gene-barcode matrix for all samples.
Single-cell RNA-Seq experiments are subject to batch effects. The CCA method from Seurat and the MNN method from scran are both available in the respective R packages to account for this (with different approaches). Batch effect correction needs to be performed before downstream analysis to ensure comparisons are valid.

This methods require an overlap between the samples to estimate technical errors and account for them. You need to bear in mind that neither of these methods works unless:

some of the cells detected in each sample are expected to be the same cell type
biological replicates of each sample have been performed

Once you have performed this correction, then standard Seurat workflows such as those shown in the tutorials can be performed which will identify clusters. You can use these to identify clusters specific to each genotype. Bear in mind that some overlap should be expected and if there isn't one, batch effect correction will over-correct for this.

Hi Tom, Thank you for your advice. You have mentioned that we can use a remote server. Yes I have tried that. But I don't know how to save our data in a remote server. Can you tell me how to? Thank you so much!

CCA corrects based on shared variability. ie. if you expect to find a unique new cell type or cell state for your mutant condition than after CCA that cell type will more often than not be blended in with the other clusters. You won't detect it. If this is what you expect MNN might be better. However if you expect changes in gene expression between the same cluster for wt and mutant CCA might work just fine. Those limitations have been addressed very recently by many groups.

MNN has similar limitations. It’s stated as such in the original publication of this method.

You just need to run cellranger, it will write output files to disk. You save data to disk from R on a remote server. This is really a separate question and I think you’ll find answers for how to do it elsewhere. You need to ask specific questions in SE to get the best answers.

There is no single ended answer to your question as everything will depend on your data. For example: Is this unsupervised data where you have to browse to 100 cell types for each sample or did you FACS sort on specific markers and are you looking for specific differences between 2 specific cell types etc...
Below you can find 4 recent packages that assist in your question. I strongly suggest you read the preprint of each, play around with your data and then follow the suggestion in the comments and ask more specific questions about how to proceed correctly.
If you want to use Seurat you can start here:
Seurat 3.0 pre-release to integrate your data sets:
https://satijalab.org/seurat/pancreas_integration_label_transfer.html
Seurat 2 tutorial on CCA to compare different conditions
you can use these methods on your integrated data sets.
https://satijalab.org/seurat/immune_alignment.html
LIGER Package and preprint can be found here
https://macoskolab.github.io/liger/
HARMONY
https://www.biorxiv.org/content/early/2018/11/05/461954
Conos
https://github.com/hms-dbmi/conos