text:
programmatic secondary structure prediction for >36-mer DNA oligonucleotides
I'm writing a tool to automate Sanger sequencing primer design for a production lab that uses a universal-tail chemistry Sanger sequencing to verify NGS results. Essentially, the template DNA is amplification using primers with universal 5' tails (M13F/R), and then the "BigDye" reaction is performed with an M13F and M13R primers in separate wells.
I'm able to design logical primers for each region in a request using Primer3 (v2.5.0), but then I'd like to calculate folding thermodynamics for the M13-appended primers to determine the best primer partner for the 5' M13 appendages.
In other words, in the last step, I would like to know if I should use:

M13F-forward_primer + M13R-reverse_primer

or

M13R-forward_primer + M13F-reverse_primer

The problem is that Primer3 won't fold an oligo > 36 nt, since the algorithm for Tm is diminishingly accurate for longer oligos. However, I believe thermodynamic predictions are valuable, even if only relatively, i.e. comparing one 5' appendage to another 5' appendage on the same allele-specific sequence, so I'd like the most painless answer, even if it isn't accurate in absolute terms.
I'm aware that UNAFold will do this, but it seems like extreme overkill for what I need to do.
For lack of a better solution, I've implemented a workaround wherein I pair each allele-specific primer with either M13F or M13R and submit the pair as a check_primer task. In other words, the left primer is one of the allele-specific primers that primer3 designed, and the right primer is either M13F, or M13R. I then choose the M13 appendage orientation based on the PRIMER_PAIR_0_COMPL_ANY_TH output values for the two possible orientations.
I would like to know: Is there a better solution? Primer3 docs say that the Tm calculation (and not necessarily folding thermodynamic calculations) is inaccurate for 37+ mers. Is there a way to force Primer3 to fold a sequence > 36 nt? If not, is there a simple tool for folding >36-mers?
Follow-Up: The workaround described in my OP was inadequate and UNAFold would cost $9k (which would take months for my employer to approve - don't get me started), so I tried @MichaelG's solution below. It turned out to work quite well. I'm computing the maximum Tm value for all fragments from every oligo, and that number agrees quite well with the browser-based UNAFold output.

Welcome to the site. Do you want to have a script or another program to calculate it or force somehow Primer3 to calculate it? I doubt we can force Primer3 to do this calculation but it would be easier to help you with your script for thermodynamic predictions... If you could clarify why UNAFold is an overkill it will be helpful too (from my point of view if it does what you want and it easy and reliable I would use it).

I agree with @llrs Primer3 is great but if the error limit of the algorithm has been reached there is no point pushing it any further. If you really wanted to down that path you could just obtain the source code and remove the restriction.

@llrs - I'm looking for a solution I can implement in Python or Bash - I only know those, and I'm only really adept at Python. TBH, my question is partially to confirm that primer3 REALLY won't do a >36-mer, because it seems like a significant limitation for an otherwise great tool. If that's the case, you're right that UNAFold isn't overkill. See my description of a workaround in Primer3 - I think that's the route I'll take for now - I'll implement UNAFold if we can't get sequence data from our test set.

@MichaelG. I understand your point, but I just want to clarify: the p3 documentation says Tm calculations are inaccurate for long oligos, but I want the results from folding thermodynamic calculations.

If you really insisted on using Primer3 you could perform a 36bp sliding window analysis. You would submit your primers in 36bp 'chunks' with a step size of say 3 bp. You would need to "circularise" your primer within the sliding window to assess the "beginning" and the "end" and remove "artifactual" hairpins. Its doable via a script and you would automate the web upload and data extraction via a "digger/spider" or else just use a download equivalent (I've used a standalone). There would be a "blindspot" using this approach, which you could overcome if you wanted.
So summary, its complicated and I think your suggestion of using UNAFold is a better idea.

I'm not insisting on primer3, but that's a really interesting solution. I came up with a different workaround after I posted this question - now incorporated in an edit to in my original. I'm going to test that against the answer I would have gotten from UNAFold (manually via UNAFold's browser interface) to see if I would have made the same choice with UNAFold data. I think that will tell me whether I need to invest the time into implementing UNAFold.


text:
How can I remove (non-trivial) duplicates from a VCF file?
This is related to the question I asked here. Consider a vcf file that contains duplicate variants, but where the duplicates aren't simply the same thing in the same notation but instead one is a subset of the other. For example:
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

These the two variants are actually the same. They result in exactly the same genotype. Changing  AACAC to AATAC at position 529514 just means change C to T at position 529516.
Is there any tool that can detect such duplicates and remove them? I tried vcfuniq from vcflib, but that doesn't seem to recognize this as a duplicate. I think it only looks at the 1st 4 fields and only considers duplicates those variants with exactly the same values in the 1st 4 fields:
$ ./bin/vcfuniq test.vcf
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

However, as explained in the linked question, EBI's vcf_validator considers this invalid. And it doesn't really make sense to have these duplicates in any case, so is there any way I can detect and remove them? Preferably an existing tool, but I am open to scripting solutions as well.

This is further complicated by cases like this one:
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC       0.00    .       .     GT    0/1
chr12   529516  529516  C       T       .       PASS    .     GT    0/1

Unfortunately, this one won't be caught by the approach in Daniel's clever script:
$ cat test2.vcf | foo.py
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC   0.00    .   .   GT  0/1
chr12   529516  529516  C   T   .   PASS    .   GT  0/1

It turns out that bcftools can do this (tested on bcftools-1.8), if you give it the reference genome to test against:
$ bcftools norm -d none -f hg19.fa test.vcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -d none -f hg19.fa test.vcf; Date=Wed Feb 27 16:08:44 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  .       C       T       .       PASS    .       GT      0/1
Lines   total/split/realigned/skipped:  2/0/1/0

For the more complex case of the multi-allelic variant in the second VCF example from the question, you can run it through bcftools twice. Once using norm to left-align and split multi-allelic variants, and then again to remove the duplicates:
$ bcftools norm -m -any -NO z -O v -o - ~/test2.vcf |
    bcftools norm -d none -f hg19.fa
Lines   total/split/realigned/skipped:  2/1/0/0
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -m -any -NO z -O v -o - test2.vcf; Date=Wed Feb 27 18:18:32 2019
##bcftools_normCommand=norm -d none -f hg19.fa -; Date=Wed Feb 27 18:18:32 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  529514  CAC     T       0       .       .       GT      0/1
chr12   529516  529514  C       T       0       .       .       GT      0/0
Lines   total/split/realigned/skipped:  3/0/2/0

Great information to have in the back pocket. Nice solution!

Wish I could +1 this again.

I'm no expert with VCF (few can say they are!) but I have worked a lot with VCF data in the last few years, both tools to consume and produce VCF. I've never seen variants encoded in this fashion, and it seems to be non-canonical. Typically:

Single nucleotide variants (SNVs) are encoded with a single base as the REF allele and a single base as the ALT allele.
In the case of insertions or deletions, the shorter of the REF and ALT alleles will be a single base, the base preceding the inserted/deleted sequence. Thus the first base of the REF and ALT alleles is always the same.
In the rarer case of two or more consecutive substitutions forming a multinucleotide variant (MNV) the REF and ALT alleles will have the same length.

Using multi-bp strings of the same length to encode SNVs is unnecessary and, as you've pointed out, problematic. This makes me think its a bug or a "feature" of the variant predictor that produced the VCF.
In this case, I'd write a small script that would check for variants where the REF and ALT alleles have the same length. If the base is the same for REF and ALT in any position, drop it, and adjust the position accordingly.
The script below will convert these funky SNVs to the canonical representation, and will also work on MNVs. Standard tools should then work to remove the duplicates.
#!/usr/bin/env python3

def canonicalize(instream):
    for line in instream:
        if not line.startswith('#'):
            values = line.split('\t')
            pos = int(values[1])
            ref, alt = values[3:5]
            if len(ref) > 1 and len(ref) == len(alt):
                # How many bp to trim off the end
                for n, (r, a) in enumerate(zip(ref[::-1], alt[::-1])):
                    if r != a:
                        revoffset = -1 * n
                        break

                # How many bp to trim off the front
                for n, (r, a) in enumerate(zip(ref, alt)):
                    if r != a:
                        offset = n
                        values[1] = str(pos + offset)
                        values[3] = ref[offset:revoffset]
                        values[4] = alt[offset:revoffset]
                        break
                line = '\t'.join(values)
        yield line

if __name__ == '__main__':
    import sys
    for line in canonicalize(sys.stdin):
        print(line, end='')

UPDATE: Upon further reflection, the more complicated example you listed actually makes sense. In the reference we have the sequence AACAC, and the alternate alleles represent two variations on this: deletion of the final two bp (in the first case), and a point mutation of the middle C to T (in both cases). Usually, only a single bp precedes the definition of an indel, so I would've encoded this complex variant as ref=ACAC alt=AT,ATAC.
So the SNV is "implied by"/"encoded in"/"redundant with" the complex variant, but it's not strictly a duplicate. I'm curious whether the VCF validator complains about these cases as well?

It is indeed atypical, but this question was prompted because I actually encountered this in the wild. I didn't generate the vcf, it was sent to me, but the header suggests it was produced by freebayes and the merge of two separate files. So this was probably an artefact of the merging. Unfortunately, I need to deal with VCF files that are given to me by clients, so while I can insist that they conform to the standards, this does conform (AFAIK), so I needed a way of fixing it.

Nice one! Unfortunately (and sorry, I should have made this clear) the actual file I had contains multi-allelic variants where only one was a dupe (see updated question). Your approach won't catch those, but bcftools does (see my answer). This works great for singles though.

Dude, that is some gnarly VCF.

Welcome to my world :(

Still a cool solution @DanielStandage!

source: Stack Exchange

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

sort: votes


text:
scRNASeq expression matrix with decimal values
I am trying to replicate some results of a scRNASeq experiment and, when I looked at the data provided by the author, I noticed that some of the counts in the expression matrix are represented as decimals.
Would this be because the authors uploaded a normalised count matrix or for other reasons? Also, please note that the authors used the 10X genomics cellranger pipeline for the processing of raw data.
I have provided a brief snippet of the data frame here:
structure(list(P1TLH_AAACCTGAGCAGCCTC_1 = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), P1TLH_AAACCTGTCCTCATTA_1 = c(0,
0, 0.314759545319035, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.314759545319035, 0, 0), P1TLH_AAACCTGTCTAAGCCA_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18024827050334, 0, 0, 0, 0, 0, 1.18024827050334,
0, 0, 0), P1TLH_AAACGGGAGTAGGCCA_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.644905372424165, 0, 0, 0), P1TLH_AAACGGGGTTCGGGCT_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.780522072779651,
0, 0, 0), P1TLH_AAAGCAACAGTAAGAT_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), row.names = c("RP11-34P13.7",
"FO538757.2", "AP006222.2", "RP4-669L17.10", "RP5-857K21.4",
"RP11-206L10.9", "LINC00115", "FAM41C", "RP11-54O7.1", "RP11-54O7.3",
"SAMD11", "NOC2L", "KLHL17", "PLEKHN1", "RP11-54O7.17", "HES4",
"ISG15", "AGRN", "C1orf159", "RP11-465B22.8"), class = "data.frame")

Also, I am quite confused of the structure they used - as I don't have much experience with single cell datasets. Just to give some further context:
The colnames of the "full" dataframe are coded from P1TLH_sequence_1 to P5TLH_sequence_1. There are 5 patients in the study.
I'm assuming that each of the colnames represents a cell for a particular patient while the rows are the genes.
If this assumption is correct, how would I do the conversion from colnames in the dataframe to the "cell names"? And lastly, why are there decimals in this expression matrix?
Thank you all for your time,

It looks like you have normalized counts. My guess is that they are log-normalized with some scale factor, probably 10000.
With the whole matrix you should be able to work out the scale factor for each cell. See https://www.nxn.se/valent/2018/10/25/unscaling-scaled-counts-in-scrna-seq-data
They typical format for scRNA-seq data is to have cells in the columns and genes in the rows. So here your rownames are the gene names and colnames are the cell names.

Thanks for the response ! Do you have any idea of a DB of some sorts that is used to convert the colnames from sequence format to cellID format?

What's the cellID format?

It is not a particular format called cellID per ce. I just mean conversion from the sequences to a naming system that would be more suitable for downstream analysis / clustering.

Ok in that case any unique identifier is fine and I would leave the cell names as they are. The only thing you might like to do is extract the patient ID information that appears to be encoded in the cell name. If you use Seurat, this will be done for you when creating the Seurat object.

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


text:
Plink 1.9 --merge-list changing order of fam file
I am doing some whole-genome sequence analysis and to speed things up split the data into chromosomes and also into regions within those chromosomes. (So I have my dataset split into ~300ish files) I am doing some IBD analysis with plink files and to do this IBD analysis I am combining the files into one big plink file.
I ran into issues today where when using plink 1.9 and the --bfile and --merge-list commands, the output fam file is different from the input fam file.  My dataset has related individuals so the fam file is very important.  Also its is important to note that when using binary files the order in the fam file needs to be exactly the same as the order the individuals appear in the other files (I learned this the hard way haha).
Does anyone know why plink 1.9 changes the order in the fam file when merging plink files, and how to fix this?  It works with plink 1.07 but with 1.9 I get a different fam file.
Performing "diff" on input and output fam files shows they are different:
plink=$mydir/bin/plink_1_9/plink

$plink --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

Using the same code but using plink-1.07 creates matching fam files as expected"
plink1=bin/plink-1.07-x86_64/plink
$plink1 --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

The merge_list.txt file is in this format:
$indir/chr22/286/$prefix.chr.22.region.286.qc.bed $indir/chr22/286/$prefix.chr.22.region.286.qc.bim $indir/chr22/286/$prefix.chr.22.region.286.qc.fam
$indir/chr22/287/$prefix.chr.22.region.287.qc.bed $indir/chr22/287/$prefix.chr.22.region.287.qc.bim $indir/chr22/287/$prefix.chr.22.region.287.qc.fam
$indir/chr22/288/$prefix.chr.22.region.288.qc.bed $indir/chr22/288/$prefix.chr.22.region.288.qc.bim $indir/chr22/288/$prefix.chr.22.region.288.qc.fam

Can someone explain to me why this could be happening and if there is a way to fix it? Thanks so much!!
(Note, I know I can manually replace the fam file with the correct one with cp or mv, but I would like to know how to handle this with plink 1.9, if possible, as I am afraid it could be indicative of an underlying error in the merge process :) )

Thanks Christopher Chang for the great answer via the plink2 google group! (See here)
His answer was as follows:
"Yes, plink 1.9 did change the default merged sample order.  However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge.  (You can also use --indiv-sort + --make-bed at any time to change the sample order to whatever you want.)
As for why the default was changed, take a look at the new default order: it will probably make a lot of sense.  It’s “natural-sorted” which tends to make it easy for humans to manually search through regardless of what mix of letters and numbers are used in the IDs."
See also the Plink 1.9 documentation for --indiv-sort which makes the point that this behavior is default only for merges--in other commands Plink 1.9 uses the same ordering as in Plink 1.07, which is the order the samples were in.


text:
Why does my SplitDotPlotGG plot show differential expression of a gene but Seurat::FindMarkers() function does not return it?
I have analysed  and clustered my single cell rna seq data with methods in the Seurat package. After clustering, I obtained two main clusters: 0 and 1. I used the Seurat SplitDotPlotGG function to visualise the differential expression of a gene of interest X.
SplitDotPlotGG(object,
               grouping.var = "clusterID",
               genes.plot = c("geneX"),
               cols.use = c("blue", "red"), plot.legend = T)

The plot shows clearly that there could be a significant differential expression of geneX between the conditions represented in clusters 0 and 1. However, I cannot see geneX in the list of differentially expressed genes returned when I run:
FindMarkers(object, ident.1 = 0, ident.2 = 1, genes.use = "geneX")

Is there be any explanation for this?
Thanks in advance.

You're showing groups 1 and 2 in the plot but comparing groups 0 and 1 in the differential expression test, so you need to instead set ident.1 = 1, ident.2 = 2 in FindMarkers.
The FindMarkers function by default only returns genes with >0.25 log fold-change expression, with >3 cells in each group. It's possible that these criteria are not being met for gene X. You can try changing these defaults, setting logfc.threshold = 0 to return all genes.


text:
Demultiplex nanopore reads with custom barcodes
We have a problem trying to demultiplex MinION sequences with custom barcodes. Do you have any software recommendations we can try for demultiplexing or how to demultiplex these custom barcodes with Albacore? We have tried using albacore but it only recognizes barcodes provided by Nanopore Technologies.

I've written some Perl code for demultiplexing nanopore reads using LAST. A workflow demonstrating the use of this script (and how it works under the hood) can be found here:
https://dx.doi.org/10.17504/protocols.io.14egnxw4zl5d
You feed the script a barcode file and input fastq file, and it creates a trained mismatch matrix, aligns to the barcodes, then demultiplexes the reads into separate files for each barcode:
./fastq-dental.pl -barcode barcode_full_PBK004.fa -mat bc.mat reads_all.fastq.gz

Example output:
         9 BC01
        39 BC02
         1 BC02_BC07
        16 BC03
        38 BC04
         2 BC04_BC06
         1 BC04_BC07
         1 BC04_BC09
       643 BC05
         2 BC05_BC07
    493175 BC06
       115 BC06_BC06
         2 BC06_BC07
    663161 BC07
         1 BC07_BC05
       117 BC07_BC07
         1 BC07_BC09
       106 BC08
         1 BC08_BC06
        75 BC09
        72 BC10
         1 BC10_BC07
        41 BC11
        44 BC12
         1 BC12_BC06
     13535 BCchim
    237796 BCnoadapt
       152 BCnone
        90 RB12A

While the approach designed around using the ONT barcodes, the barcode fasta file can be replaced by another file [using barcode ID in the header] without change the protocol.
FWIW, there is also now some support for detecting strand direction (e.g. using the strand switch primer provided by ONT), but I haven't yet got it working properly doing both demultiplexing and strand orientation at the same time.

Thank you! I am using minibar, works great and very easy to add custom barcodes.
L.Y.

There are two other ways in which you could demultiplex reads if guppy doesn't have the config file for the kit you are using.

Minibar: This works for dual index barcodes specifically
Porechop: This is intended for the ONT sets, but you can also replace them with your own. However, it is designed to work with single index barcodes. That means, you will have to run it once with your forward indices, then take that output and run again with the reverse, and then do some gymnastics to combine it all together

What do you mean with run it twice, with forward and reverse indices?


text:
Normalization for two bulk RNA-Seq samples to enable reliable fold-change estimation between genes
I have two bulk RNA-Seq samples, already tpm-normalized.
I would like to know what is a reasonable normalization procedure to enable downstream log fold-change estimation.
The distribution of the two samples using the common set of genes looks similar:

However, the two samples have only been tpm-normalized, is it enough to guarantee reliable fold-change estimation? Should I use another normalization procedure, e.g. Quantile Normalization, before comparison?
My objective is to define a signature using the genes that are up-regulated in Sample1 with respect to Sample0, and vice versa. I'm using log fold-changes, but I'm concerned that their value may be affected by each sample distribution.
Do you also have suggestions for the definition of up-regulated genes with these data?

You have only two samples?
You aren't going to be able to draw strong conclusions from that no matter what you do.  Clever statistics don't work without replicates.

It's not a good idea to do tpm normalisation prior to differential expression analysis, because the actual read counts are useful to determine shot noise and statistical significance. DESeq2 includes read normalisation as part of its methods for differential expression analysis.
I think shot noise is best explained in terms of shooting photons at a target, as found on Wikipedia:

Shot noise exists because phenomena such as light and electric current
  consist of the movement of discrete (also called "quantized")
  'packets'. Consider light—a stream of discrete photons—coming out of a
  laser pointer and hitting a wall to create a visible spot. The
  fundamental physical processes that govern light emission are such
  that these photons are emitted from the laser at random times; but the
  many billions of photons needed to create a spot are so many that the
  brightness, the number of photons per unit of time, varies only
  infinitesimally with time. However, if the laser brightness is reduced
  until only a handful of photons hit the wall every second, the
  relative fluctuations in number of photons, i.e., brightness, will be
  significant, just as when tossing a coin a few times. These
  fluctuations are shot noise.

Sequencing also has shot noise effects because the sequencing process is a random sampling method. For low-abundance targets, the likelihood of sampling a target is low, so the effect of a single hit is amplified. This reduces the statistical significance when comparing different low-abundance expression values, because the variation in expression is highly dependent on which targets won the sampling lottery. The shot noise variation / expression curve can be seen in Figure 7 of this paper.

I agree with TPM for a lot of reasons, unfortunately the data was already in TPM. Can you explain more about how read counts are useful to determine shot noise and statistical significance? Thanks

What I have generally done in the past is to process the data using voom in the limma package for bulk RNASeq. Inside voom you can call for different normalization methods to be used - "TMM" works fine for me and, is advocated by many in the field.
voom will output an object containing the normalized expression values in a log2 scale, which, also in my experience, has worked out just fine for calculating log fold changes.
Check out this link for more info on the package as well as normalization methods: https://www.bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html
It is a very thorough introduction to the package and all of its capabilities.
Good luck!


text:
Why are the clusters in the left side of a tSNE plot and not in the centre?
I have no idea why I used the subsetdata in Seurat after using CCA the TSNEPlot shows the cluster all in the left side and merge， did I miss something? Here are the code and figure I got.
cluster <-SubsetData(object = immune.combined, ident.use= "3")

cluster <- NormalizeData(object = cluster, normalization.method = "logNormalize", scale.factor = 10000)

cluster <- FindVariableGenes(object = cluster, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0, x.high.cutoff = 6, y.cutoff = 0, do.plot = F)

cluster <- ScaleData(object = cluster, display.progress = FALSE)
cluster <- RunPCA(object = cluster, pc.genes = cluster@var.genes, pcs.print = 1:10, genes.print =5)

PCAPlot(object = cluster, group.by= "orig.ident", dim.1 = 1, dim.2 = 2)

cluster <- FindClusters (object = cluster, reduction.type = "pca", dims.use = 1:10, resolution =0.3, print.output = 0, save.SNN = TRUE, force.recalc = T)

TSNEPlot(object = cluster, do.label = TRUE,do.return = TRUE,pt.size = 2,label.size = 5, no.legend=FALSE)

I'm not familiar with Seurat, what exactly is the data your are performing your tSNE analysis on? 50 principle components of your differentially expressed genes of scRNA data?

The lack of data transformations in Seurat is worrying ("logNormalize" is the only one). The idea of "one size fits all" is never cool in biology. I don't know Seurat RNA-seq, but I know the surrounding statistics. Have you done a PCA plot?

I mean, we can't really tell how dense those isolated clusters are in the other three quadrants. How many cells do you have? tSNE is stochastic and depends heavily on parameters, so default settings aren't good for every data set. I'd recommend having a look at this blog post, and then rerunning using several different parameter settings.

Sometimes low quality or outlier cells can "push" the rest of your data to one side of a tSNE plot. You would probably see a pattern like this in one of the principle components as well.  I would try to isolate the cells and look at their expression and if they do look like they are low quality you can likely remove them from your analysis.

I don't see CCA anywhere in your code. Is that a typo? Should be PCA?

The issue here is simply that you did not re-compute the tSNE after subsetting the data, so you're looking at a subset of cells from a larger embedding. You would never get an embedding that looks like that by running tSNE on only those cells.
You also don't need to re-normalize the data after subsetting the cells. Since this is a per-cell calculation, the results will be identical to those in the larger matrix (you can check this for yourself).

Hi TimStuart, I have tried no normalization and no scale and found out it is no a difference between the one I scale and normalization them. Do you know what's the problem it is?@TimStuart

Scale data will be different, normalize data will not be different. You need to compute a new tSNE: cluster <- RunTSNE(object = cluster, reduction.type = 'pca', dims.use = 1:10)

I'm going to stick my neck out here. Overall PCA followed by tSNE is a very powerful method of unsupervised learning however, I agree with @GWW,
cluster <- ScaleData(object = cluster, display.progress = FALSE)

If this is a max-min rescale transformation (common in supervised learning) then the outliers will badly affect the clustering. You can remove outliers that are over 2 sd from a formal analysis and re-run the calculation. That is just the risk of a max-min rescaling.
The one issue I'm concerned about is the initial quasi-standardisation followed by log-transformation,
log1p(value/colSums[cell-idx] *scale_factor)

In this case you are multiplying by 10000, this is to remove the ratio (used as a quasi-standardisation) and then taking the log ??? That is not the way to transform a ratio, it would be better to take the reciprocal of the ratio. There are other (better) methods, a formal standardisation would be useful, before a transformation is considered. Possibly if that was done correctly you wouldn't have the problem of outliers potentially crunching your data.
If you fiddle with the parameters you might get something - but its a symptom not the cause of the problem. However, if removing the outliers works (I have not checked the ScaleData function), then you are sorted and all else is forgotten (at least for this data set). The information I read was here.

What is the issue with log1p(counts_per_10k)? This is a common way to normalize for different sequencing depth between cells

If the scale factor is based on an underlying biological reading (cell count) -- no problem at all. To multiply by 10000 default suggests this isn't a standardisation (because the scale factor would be variable). Standardisation against a control is a better approach IMO to offset experimental variation between runs.

The scale factor is variable, it is the total counts for the cell (ie colSums(cell)). We then multiply this by some large constant value (1e4,5,6 etc) to avoid dealing with small numbers

It doesn't appear to being used like that in the above code. The default value set in the documentation was the default value used by the OP

Sorry this is a bit confusing because I used the term "scale factor" above referring to the cell-specific scaling factor, not the argument in the function. The argument scale_factor in the function refers to the constant scale factor (by default 10000). A cell-specific scale factor does not need to be specified as it is calculated as the sum of counts in the cell. In the code above, colSums[cell-idx] is the cell-specific factor and scale_factor is an arbitrary constant set by default to 10000


text:
Why does a missing label 11plex TMT shows up at almost 50% intensity compared to other labels?
We have a test plate that was done with 10 of the 11plex TMT labels. So, the intensities of the 11th label should be almost zero, plus the leakage from the 10th label due to 13C content. Although we observe a substantial drop in the intensities, we see ratios up to 50 percent compared to the highest reference intensity, which is the maximum of the intensities of the other labels, if we sum up the intensities for all proteins. For individual proteins the intensity sometimes is higher than for all "real" samples. Only talking about proteins that are identified "Confidential". Is that something you would expect? Is this due to the normalization?
I was using Reporter 0.7.17 from command line without specifying a reference sample. Pretty much default configuration and based on the default protein report.

Fig: Ratios (intensities) reported by the program.

Wow, thats a tough problem. Perhaps use a reference and standardised against it (???). You might then get something.

The developers answered along these lines:
This might be a normalization artefact. Reporter is normalizing the ratios so that the median is 1. So, if the label is almost not present the noise gets increased. To check that you can turn of the normalization. This function is currently not documented. The developers are working on it though.
GitHub issue


text:
GBS: clustering of PCoA axes with R package mclust to describe identity by missingness
I have a SNP dataset with 56 populations and 430 samples that was produced by GBS sequencing. I need to identify groups of samples that have similar patterns of genotype missingness, since there is substantial structure in missing genotypes, shown here in a PCoA ordination of genotype matrix with 1-0 values, where 0 represents missing gt and 1 is presence of a gt value:

Here the missingess threshold for genotype inclusion is <20%.
I want to cluster the samples using the scores on 20-30 PCoA axes, which is 17%-21% of variation, using the R package mclust. I ran ape::PCoA() with an Euclidian distance matrix from paralleldist::pardist(). When I run mclust::mclustBIC(), I get a weird result in that the there is no modal value for any of the better-supported models, here shown in the legend:

I had expected to find a model that would support the existence of a small number of clusters (2-6 maybe).  I also don't understand why the VVE model only has values to 3 clusters, where it presents the maximum value. Any ideas?  I have tried relaxing the minimum missing level, and using data with differing minimum locus depth levels.  I get very similar results.
Could someone please help me understand what is going on here and how I might address it?  Once I have clusters, I can look at the assignment of samples to clusters by population, flowcell, and other stratifying variables. I would be happy to edit to include more information, if that would be helpful.

So, what do you want help with: Why there isn't a single model with a definite number of models? Why VVE model has only values for 3 components? or How to identify subpopulations? How many SNPS do you have? BTW it would help to visualize better the PCoAs if you use scale = "free" in facet_grid

Hey @llrs following our earlier discussion ML could actually work here, like a random forest. Here, the model assumptions don't conform to the data, so rather than figure out the "magically distribution" you just wrap ML around it. There is a caveat (PCA isn't a nice ML approach in this case because it would be hard to interpret the result), but it would 'solve' it.

@MichaelG. I'm not sure that model assumptions don't fit on this data, but that the result is not the one expected which is a completely different thing in my opinion. I don't understand why there should be several subpopulation in first place. The random forest has its own assumptions that should be checked too,  I don't understand how using it here would help here.

@llhs that is certainly a possibility, I did think about it, but either way the BIC must reach a plateau. Random forest, if used with a PCA, doesn't mean very much. If you have an a priori cluster and you want to know (in this case) what SNPs contribute to the sub-population, a random forest is one approach that will work (not using PCA). This analysis per se is interesting and if it worked, could useful. Honestly I would I simply use a traditional approach in pop gen first. PCA is when stuff doesn't work, but only in pop gen.

The final thing to add is "missing genotypes,", if this means that the experiment didn't work for that particular allele, its no biggy. If it means there was an expected genotype and another genotype was found in its place, then that is a big deal.


text:
Understanding PCHeatmap outputs
I am currently trying to understand the purpose of these PCHeatmaps - part of the seurat package in R:

All the online documentation I have searched for has only highlighted how it is used to observe heterogeneity in the data.. This much I understand.
My questions are:

Are the genes in the rows of each PC plot just the most variable x
number of genes defined in FindVariableGenes?
If not, then how are each of these genes 'chosen' for use in the PCHeatmap
visualization?
What would be the different conclusions/deductions drawn from the PC1 heatmap and PC2 heatmap with regards to heterogeneity for example?
And finally, how can these conclusions help in the analysis of scRNASeq data?

The PCHeatmap function (renamed DimHeatmap in Seurat v3) can be used to help determine the number of principal components to use in downstream analysis, as well as to visualize the top genes contributing to each PC.
Both cells and genes are ordered by their PC scores, and by default the 15 genes with the highest and 15 genes with the lowest PC loadings are displayed for each PC.
In your case with PC1 and PC2, my conclusion would be that PC1 and PC2 both strongly separate different populations of cells. For PC1 there is a smaller population with very strong separation, and PC2 splits the cells roughly in half.
You can also set the cells.use parameter (renamed cells in Seurat v3) to only show the top n cells for each PC which can help with visualization.

What do you mean by PC score? Do you do two PCAs where one uses the genes as features and the other the cells? I don't get it.

@Inhibitor no just one PCA. but roughly speaking the cell scores are the eigenvector values, and the gene score is the linear combination of the cell score and the cell counts

@august Thanks. I don't have the intuition to immediately see how these things become useful scores (need to think about it a bit), but at least it's something I can calculate. Wish people would write this in their articles, coz' if it confuses me then it certainly also confuses wet lab biologists...


text:
Is nanopolish worth it since faster polishing software is available?
For Oxford Nanopore contigs produced by any long-read assembler has anyone performed any benchmarks to compare the polishing tools racon, nanopolish, or something like wtpoa-cns?
The compute time of nanopolish is quite long and I would rather use another tool if it will suffice. After long-read polishing I will polish more with two rounds of pilon using 200x coverage Illumina shotgun data from the same individual.

Nanopolish is necessary if you want to get high-quality consensus. racon etc don't use signal data. They can't achieve high quality. pilon at times doesn't work well. Even if it worked perfectly, it wouldn't help repetitive regions inaccessible by short reads.

You're quite right. I polished twice with racon then twice with pilon. Mapped short read pairs and long reads. All of the gaps in the short read mappings were covered with long reads. Going to use nanopolish now!


text:
Repetition-free Longest Common Subsequence as similarity measure
I'm starting a thesis in Computer Science with topic the RFLCS (Repetition-free Longest Common Subsequence). I know very little about Bioinformatics problems and Genomics. I would like to understand the motivations behind this problem.
More precisely, I'd like to understand why one need a similarity measure between genomes, and how the RFLCS is an useful measure.
[Adi et al.] introduce the problem in this article.

What is the RFLCS? PLease [edit] your question and give us some more detail so we can understand exactly what you're asking. However, I am afraid that this will probably be too broad to answer (see the "too broad" section here: https://bioinformatics.stackexchange.com/help/closed-questions).

Done! I narrowed the field, hope enough.

Still a bit obscure. I think the question is there to sharpen your coding skills. Its fairly easy to solve. Although it could be solved as a 5% tail of a Poisson distribution (I doubt that is what its about).

I found the paper introducing it and I think you should mention it in the question. I think that your supervisor is the best person to talk to regarding motivations of the problem. If there is anything you would like to get help with from the community, you need to be more concrete. "Why do I do my thesis" is way too broad.

A brief literature review can give you the answers.
This preprint gives the definition of RFLCS more clearly:

A repetition free Longest Common Subsequence (LCS) of two sequences x and y is an LCS of x and y where each symbol may appear at most once.

Adi et al proposed to solve the genome rearrangement problem raised in this paper. If you want to know more, follow the references and citations of these papers and read them.
Basically, Adi et al think the length of RFLCS measures the distance between two genomes. The longer the more similar. I am not sure if this is a good measurement, though.

Now I see. If shifts in LCS were solely the consequence of rearrangements .. yes, but we all know thats not true. Transposons move, microsatellites, minisatellites and DNA virus integration (e.g. HPV) will all be confounders,

Thanks for the answer. I read several papers, but no one explains why RFLCS should be a good measure. Seems nobody have assessed its effectivness.

@DanieleCuomo I don't think RFLCS will be effective, tbh.


text:
I have 23andme text files and would like to convert to SAM/BAM format
I would like to convert 23andme text file to NextGen to BAM file for Yfull.com to read.  It is difficult to get answers on how to convert to SAM/BAM file for the 23andme text file to be converted to SAM/BAM tools

Please include an example of the text file since 23andme is not a standard format or something that people in the field will necessarily be familiar with. But if you mean you just have a list of variant positions (which is what this sort of consumer service usually provides), it is impossible to convert to bam.

I would be happy to provide a text of 23andme. Is there a way to upload it?

We just need to know what the file looks like so we can understand what you need. Just copy the first few lines (if there is a header, make sure you include both the header and a few lines of data) and paste them into your question. Use the formatting tools to format them as code.

It would also be helpful if you could explain what your final objective is. BAM files are used to store the alignment of raw sequencing data (reads) against (usually) a reference genome. 23&me should have given you a set of genomic variants, something that is produced by analyzing a bam file. Why do you want to go back to a bam file now?

I want to upload the Y-DNA to Yfull.com to NetGen for any markers 23andme didn't test all markers; Yfull.com does not read text files.  They read BAM files.

Please show us a few lines of your 23&me file. What you are asking for is almost certainly impossible unless 23&me have given you the raw reads. We can't know unless you show us a few lines so we know what you have.

Here is the data for 23andme.  I am not sure how to use formatting tools. I tried to copy and paste the text file.

Note you can just [edit] your question and add/change whatever you need. You probably want to use "code" mode for the sample of your file, which you can achieve by copy-pasting the text, highlighting it with the mouse cursor and clicking on the "code" formatting button.

Having done 23andme myself I can tell you that your variant file, which contains SNP genotypes, cannot be converted to a bam file. It does not contain the same information as a bam file. It may be helpful to familiarize yourself with those filetypes and the technologies used to obtain them. A SAM/BAM file contains alignments of reads obtained by sequencing. Your standard* 23andme analysis doesn't have reads, as it uses a SNP array for determining the alleles at predefined positions. Markers which are not included in your 23andme file cannot be added by other means because those were not genotyped by the array. Therefore also Y-DNA cannot use your 23andme file or add additional markers.
*For the sake of completeness: I am aware that under rare circumstances 23andme can perform WGS for additional variant typing, but I do not assume this is the case for OP.

If you just want to get the markers which were not reported/genotyped by 23andMe (as you mentioned in the comment), you don´t want to convert your 23andMe data to SAM/BAM (which is most probably really difficult or impossible, as explained by @Wouter De Coster), but rather to perform imputation on your data file. This will give you a genotype file with genotype or allele probabilities for all the markers. While it will not give you a 100% certainty for the lacking variants, it can certainly give you more information on many of them.
NB: you can imputate the Y-chromosome variants only if your 23andMe file contain data from Y-chromosome. If you want to know about that option, let me know.

I think the closest approximation for a SNP chip would be the intensity values.  When I asked 23andMe about obtaining these, the answer was that they don't plan on adding this feature.
Some programs use the .vcf file format instead of the 23andMe format, but I think the 23andMe format is compatible with relatively more 3rd party software than some other vendors.  If I remember correctly, I think DNA.land even has an intermediate .vcf file, which you can download.
I actually think it is important that 23andMe doesn't impute the variants in your genotype table.  However, if you are trying to dig deeper into some SNP chip data, I believe you can still get genotyping for free from Genes for Good.  They provide a few different "raw data" formats (such as .vcf files, with and without imputation, as well as a "23andMe format" file) and I got the impression that they would consider providing an even more raw form of data (to test the effect of generating genotypes on your own, even thought I don't think most sites will change).  However, to be clear, they also do not currently provide intensity values (or some other even more raw form of data).


text:
Getting genes specially up or down regulated
I have 6 RNA-seq samples like this
4 patients (005, 036, 121, 013)

I have 3 tumour samples and 3 cancer models (organoid)

This is PCA of log transformed data by DESeq2 (tumour Vs. organoid)

In this paper https://www.nature.com/articles/s41467-018-05190-9
They are saying

Finally, for each patient, the genes that are specifically up- or
  down regulated in each patient’s organoid (1 vs. all comparison, p ≤
  0.05 and abs (log2-FoldChange) ≤ 1) and that are among the top-50 genes with highest base mean expression were selected.

They have tumour and organoid like my case; My confusion is what they mean by "each patient’s organoid (1 vs. all comparison"? Because like them I need up or down regulated genes in each patient's organoid but I can not figure out I should compare each organoid versus what?
The ultimate goal would finding a signature of genes specific to each organoid as each organoid originates from tumour of a patient by these signatures one would be able to classify cancer to sub types

Nice results! It will be very easy to find the responsible genes with the tightness of those clusters You perform eg. a linear regression (GLM or machine learning) and look at the regression weights (genes). The results are so clean any method of classification will identify them. Just to mention I don't do cancer (never have, never will), but those results look cool.

Thanks a lot for your positive and encouraging comments. The assumption is that the expression pattern in each organoid should be very close and similar to the tumour from which this organoid originates. More close better an organoid would be a model of that tumour. But in my PCA we are seeing tumours are in one way and organoids are another way. My expectation was for example tumour 005 and organoid 005 placed close to each other. However I will need a signature of each oganoid as a model of cancer. You meant I take genes as independent variables and a specific organoid as response variable?

@MichaelG. Could you please elaborate on how you would use GLM to do this? Or could you send me the names of some of the packages you have used to do this? I'd like to try it out myself.. Thanks.

@FereshTeh, just came across this rather old post. Regarding "My expectation was for example tumour 005 and organoid 005 placed close to each other" part of your comment above, please see the last paragraph of the "RNA-seq analysis" section of the methods from the paper that you shared. They kick out genes that they think arise from the "artificial environment" (their own words) organoids are prone to. You would get a better looking PCA plot if you would follow their practice, however, organoid biology would still be very very (the nuance is arguable) different from tumor biology.


text:
Where can I find .cdf files for Affymetrix ChIP-chip data analysis?
I am analysing microarray ChIP data generated several years ago. I want to use the affy package for R. However, in order to load the CEL files and proceed with the analysis, I need the pombe annotation file which comes under the name of Sp20b_M_v04.cdf. I have looked everywhere for this file but I am unable to locate it. Is anyone in possession of this file or know where I can download it from?

Welcome to Bioinformatics.StackExchange @Min! Could you let us know where the instructions indicate that you need the Sp20b_M_v04.cdf file? Is this from a colleague's lab notebook, or from the affy package's documentation, or somewhere else? This info will be helpful in troublehshooting.

Hi Min, completely agree with @DanielStandage. You need to [edit] your question and add lot more details, otherwise, it's really hard to help.

@KamilSJaron I converted your answer to a comment since it wasn't really providing a full answer to the question.

@terdon I was about to delete it because it was the CEL file I linked, not cdf. ^^

@min I found the corresponding ncbi geo page. it seems that the file you are referring to isn't uploaded. You must contact authors I believe.

Crossposted with biostars?

It looks like it's in GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10187
There is apparently a zipped copy linked at the bottom in the supplementary files.
When people upload a study based on an array, there is a corresponding 'platform' record where the array details should be deposited.
Edit: you will probably need to gunzip the file first, to get it into CDF format.


text:
Unable to install bedtools on windows 10 ubuntu
I am trying to install bedtools on windows 10, but I get an error I don't understand. How can I fix it?
Building BEDTools:
=========================================================
DETECTED_VERSION = v2.27.1
CURRENT_VERSION  =
Updating version file.
 * Creating BamTools API
- Building in src/utils/FileRecordTools
  * compiling FileRecordMgr.cpp
make[1]: g++: Command not found
Makefile:36: recipe for target '../../../obj//FileRecordMgr.o' failed
make[1]: *** [../../../obj//FileRecordMgr.o] Error 127
Makefile:186: recipe for target 'src/utils/FileRecordTools' failed
make: *** [src/utils/FileRecordTools] Error 2

Hi there, I suppose you use the bash terminal available at windows 10, is that correct? The error says that you don't have a compiler for c++ called g++. You can try install it by sudo apt-get install g++. Then run make again.

Turns out that something went wrong the first time I tried to install g++, but I didn't see the error message.

My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues.
The linux version of conda installs just fine in WSL and bioconda hae pretty much any bioinformatics package you might need, already compiled.

make1: g++: Command not found

You do not have g++ installed on your machine.  This is used to compile the bedtools into an executable file.
Here is a tutorial on installing g++

If you have admin privileges, the easiest way to install g++ is sudo apt-get install -y build-essential.

@DanielStandage does that work on Windows? I assume the OP is using the newfangled Ubuntu Subsystem for Windows thing, but I don't know if that can install arbitrary packages like that.

Oops, I missed the "Windows" part and only saw Ubuntu.  It looks like apt-get is supported though. https://www.howtogeek.com/261449/how-to-install-linux-software-in-windows-10s-ubuntu-bash-shell/

I had tried that, but another error popped up without it realizing. However, the following error popped up after entering make, and I can't tell if this is critical for the functionality of bedtools:
`- Building main bedtools binary.
done.

Creating executables for old CLI.

/bin/bash: python: command not found
Makefile:145: recipe for target 'all' failed
make: *** [all] Error 127`

Please add the new issue to original post.  This is saying that you don't have python installed.  I think that if you don't want to use the legacy CLI (command line interface) you should be all set.  You can try running bedtools and see if you get a usage prompt.


text:
Significance and timing of "mux scans"
I'm using MinIONQC to do quality control on some ONT data. The software plots several read characteristics over time (hours passed during the sequencing process). These plots contain several vertical red lines. From the documentation:

Muxes, which occur every 8 hours, are shown as red dashed lines

Reading this paper (doi:10.12688/f1000research.7201.1) I was able to get a rough idea of what a "mux" or a "mux scan" is. If I understand correctly, it's a calibration step that scans each well containing 4 nanopores to determine which nanopores will be used for data collection in which order.
Do I understand the basic idea of the mux scan correctly? Does it occur at the same frequency on all ONT instruments? Does it always occur in 8hr intervals, even if the sequencing reaction is run for longer than the standard run time?

Yes, your understanding is largely correct. This originates from the situation that for each detector on a nanopore array there are 4 pores. I'll explain mux scans and groups, but this is outdated information as now another system is used.
So a MinION FC has 2048 pores with 512 sensors. Not all of these pores will be equally suitable for sequencing. At the beginning of the run originally a mux scan was performed, to test for each sensor which the best pore would be (if any) out of the 4 available. The best pores will be used for sequencing first, the rest will be turned off. That is group 1. After 8 hours the expectation is that many of these pores have degraded/accumulated damage or are not functional for other reasons. Then it is time to swith to group 2. So you change to another set of pores, which in these cases often lead to a temporary increase in yield. The default script used 8 hours, although some used adaptive software to select new pores when the first group started going bad. Also note that during the run no new mux scan was performed: the software did not check if the selection it made in the beginning still made sense. Pores in group 1 got turned off after 8 hours and switched to group 2, without taking their current 'health' into account. Cumulative yield plots get a characteristic bumpy pattern from this. Each bump is less than the previous one, eventually flatting off.
Now that's history, as the software now uses dynamic muxing and will continuously check for the best pore per channel - and as such maximize yield by having as many healthy pores as possible available for sequencing. Cumulative yield plots are now smoother.


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

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


text:
how to write the script to know each sample's cluster number?
I have followed the Seurat Tutorial:Integrating stimulated vs. control PBMC datasets to learn cell-type specific responses by using my data. And I am confused about how to write the script to calculate the cell number of each sample on each cluster.
I have run TSENPlot to know the samples and I'm using another TSENPlot to get the clusters. If I want to know find the cluster1 of each sample's number how can I write the script?

Please [edit] your question and give us some more details. What exactly do you need? Can you show us the code you have run?

Assuming you are following the basic seurat guide:
pbmc@ident

You can access the 'ident' slot of your Seurat object and you will get a vector of the cluster IDs and the names will be the respective cells.
head(seuratobj@ident)
Cell1A7 Cell1C4 Cell1D1 Cell1D2 Cell1D3 Cell1D5
0 3 3 1 3 1
Levels: 0 1 2 3 4

In my example, the printing wasn't great but essentially Cell1A7 is in cluster 0, Cell 1C4 is in cluster 3, etc...


text:
Affymetrix tags (same ID's) present in different places of the genome
could you explain me how to interpret such situation?
Same tag aligned in several distint places, near different genes...
Is the annotation able to distuinguish such 'same' tags, during DGE analysis?
Why tags arent unique?

Is the annotation able to distinguish such 'same' tags, during DGE analysis?

No, you'll end up discarding such probes (assuming you have a reason to actually use microarrays still).

Why aren't tags unique?

Platforms such as HG-U133 date back to 2002 (meaning they were designed a few years before that), when the human reference genome (and the reference transcriptome) were not nearly as "finished" as they are today. Probes that may have made sense given what we knew almost 20 years ago may no longer makes sense now given what we know. Further, microarrays are rather conservative, in that they tend to primarily gain probes over time and not have them removed (lest companies lose customers who are relying on given probes).

Thanks, for the answer! Ok, so I have the normalized expression signal from this tag. How can I distuinguish to what gene is this signal corresponds to? For example (very abstract) I have some ambiguous tag which is located in three different genes. So the signal is cumulative value of thesee 3? As a result I am not able to define the expression of these genes separately?

As I said, you'll end up excluding this probe most of the time. In modern arrays genes will have multiple probes associated with them (sometimes multiple per transcript). If you really wanted to incorporate this probe, you'd distribute it in some sort of weighted fashion to all genes/transcripts for which it's informative (in a similar to what salmon/kallisto/rsem do for RNA-seq reads). Depending on how elaborate you wanted to get that could be a lot of work, so typically you'd just exclude such probes instead.

I always assumed that these tags, such as 239963_at you show, mean that the sequence is not specific of a genome region. This means it cannot be assigned to any gene for DGE.
It could be due to changes on the reference genome or due to a faulty design of the tag sequence. However once in the microarray template they cannot be changed without introducing more problems, that's why doing RNA-seq prevent this.
Sorry I don't have any reference for this (but if someone wants to add this edit the answer).


text:
GTF from consensus sequence
I am new to using bioinformatic tools and I was hoping this community could help clear something up for me. I need to generate a gtf file. My data are a set of complete HA genes for influenza B viruses in fasta format. From reading through forums and other internet searches, my understanding is that fasta cannot be transformed to gtf or gff because these last two provide location information which fasta does not (unless available in the header). I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping.
My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence? The alignment algorithm is a simple majority rule. Because all my data are complete genes, I would think the mapping of locations (start codon and stop codon) would be the same no matter which sequence in my data set I choose to use. Is this not so?
I am looking at viral evolution in the HA segment and my end goal is to calculate the dS/dN ratio for a set influenza B viruses collected during a particular season. The program I want to use requires a gtf, but likely other options would only require the reference genome. Would it make more sense to use the strain used as the vaccine candidate as the reference genome?
I appreciate any help you can provide me!

Welcome to the site. I am not used to work with gtf, but you'll probably need a reference genome in order to align your sequences to the genome. Do you have a reference genome for the influenza B virus? What I know is that viruses change quite a lot and probably there are several genomes available. You seem familiar with alignments but you way that it is a matter of "simple majority rule", some algorithms are designed for short sequences or faster while some other are optimized for more accurate mapping, so it is not just a majority rule, but others will answer you with more detail than I can.

I think the answer is no. Let me check with my colleague who has more expertise in gtf or gff, before I cite my reasons and suggested way forward. Are you working within or associated with Derek Smith's lab (Cambridge)?

In summary, yes to this "I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping." and no to this "My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence?"

@MichaelG. Thank you for looking into this for me! No, I am not working within or associated with Derek Smith's lab in any way. Could you expand on why using the consensus sequence does not make sense?

I would still like to understand why using the consensus sequence does not make sense, but my other way forward is to use the sequence for the strain used as the vaccine candidate for that particular season. Does this make more sense? Because influenza has such a high mutation rate, it would seem that having a "standard" reference genome would not apply as it does for human genomes.

Is making a gtf file really your end goal?  I have the feeling that you are trying to do something else, and think that you need a gtf, but maybe that's not the case.

@swbarnes2 My end goal is to calculate the dS/dN ratio and the program I want to use requires a gtf. I imagine other software options would not require a gtf and instead would ask for a reference genome.

Welcome @Lindsey! It looks like these comments are getting somewhere. Please use the [edit] button above to update your original question to include a bit more about what you're trying to accomplish. It looks like you've sequenced some flu strains and would like to look at the evolution of HA genes?

I was wanting more time, but okay here's what I'm thinking,

Would it make more sense to use the strain used as the vaccine
  candidate as the reference genome?

Neither, because the GTF format is based on an HA protein structure (likely a crystal structure but could be cryoEM). This is not the same as an alignment position within a surface antigen, because most alignments of pathogen surface antigen genes will contain indels, e.g. between the GTF and your influenza epidemic. If HA doesn't have any indels fine, but I would be surprised.
The indels will cause artifactual shifts in the GTF domain boundaries unless they are corrected. Put simply an alignment position (particularly for pathogen surface antigen proteins) is rarely the same as the amino acid position of GFT because a protein structure of a surface antigen doesn't comprise indels in its native form.
If you do this work in detail, i.e. shifting between protein structure and phylogeny, you always write scripts which convert between aligned positions and "reference" (crystal/ cryoEM structure) positions. Its a common way of thinking.
So in summary, your reference sequence is the GTF sequence because that is the basis of the structure<-> alignment conversion

For positive selection analysis it does not matter what the reference is, all that matters is the correct rooting. For demarcation of domains within HA the GFT sequence needs to be at the top of the alignment 1) because it helps easily identify whether it contains indels 2) structure programs will expect the reference sequence to be there. If you are using the vaccine in a positive selection analysis ... you need to think hard about the question you are asking and the coverage of the vaccine. Epidemiologically speaking thats complicated because the passage of the vaccine will be an important source of selection and you will mistake that as the wild-type "escaping" the vaccine.

What I was suggesting was using the sequence for the HA segment of the vaccine candidate as the reference sequence used for the mapping and subsequent generation of the gtf. I'm not sure I made that clear. As you say, the HAs of the epidemic are likely to have a lot of indels, but still, some sequence must be considered the reference on which to align. In viruses with high mutation rates, what would make a good reference sequence?

I've amended above

So if I am understanding you correctly, I could use either the consensus sequence or the vaccine candidate HA sequence (or any sequence) as the reference in a positive selection analysis as long as I have the correct rooting. I'm not sure if I fully understood the last part of what you amended, but my set of HA sequences come from influenza cases, none are from the vaccine.

Hi Lindsey, .... firstly NEVER use the consensus sequence in a positive selection analysis - even if we disagree about the epidemiological focus of your study - NEVER use consequence in a phylogenetics analysis. However, yes the maximum likelihood approach to positive selection (the only one worth bothering about), the correct root is important. In other positive selection calculations (distance), it is will not make much difference.

"Transforming" Fasta files to GTF or GFF3 files is a common request, as are similar tasks such as "converting" BAM files to VCF files. But these are underdefined questions, and "transform" or "convert" isn't really the right way to think about the task.
As far as your particular request is concerned, Fasta files and GTF files contain fundamentally different types data.

Fasta files contain sequences, which could be chromosomes, assembled contigs, polypeptides, transcripts, or other biological sequences.
GTF (and GFF3 and BED) files contain intervals, indicating that some feature of interest occurs at such-and-such position on such-and-such sequence. A hypothetical example: a transcription factor binding site occurs on chromosome 2 from position 40011 to 40017.

So you are correct: it's not really helpful to think of your task in terms of "transforming" Fasta to GTF. The process of creating a GTF/GFF3/BED file is often called annotation. There are many different methods to annotate many different types of genomic features: transposable elements, segmental duplications, protein-coding genes, non-coding RNAs, and so on.
As a result, it's meaningless to say "I need a GTF file" without clearly specifying what that GTF file will annotate. (It's kind of like saying "I need an XML file" or "I need an Excel spreadsheet" without indicating what the contents of those files would be.) The things you need to clarify are:

What will the intervals in my GTF file represent?
What precise sequences will serve as the reference and provide interval coordinates for the GTF file?
How are those intervals to be determined?

It sounds like 3. is a big part of your question. If I'm reading things correctly, you want to annotate a particular class of genes on an influenza genome. Is that correct? Did you perform de novo sequencing of influenza viruses, and now hope to characterize variation in their HA genes? Stating what you hope to get out of your analysis would be helpful in finding a good way forward.

Thank you for replying! I have updated my question to be more informative of what I want to accomplish.

If I was doing this analysis, I'd want to assess selection signal against domain boundaries. In a selection calculation its one approach to demarcate the loci. For example the hydrophobic "stem" (I don't know HA terminology) will not be of huge interest in this particular question, so could be delineated to provide a control

Since you know that your sequences in your fastas are your genes of interest, why not just dummy up a gtf where every sequence is a gene?
The point of a gtf is to tell software which genomic regions matter for a certain procedure, and which regions do not.  But if you already know that your sequences are pared down to genes, then you know that the whole thing is what's relevant.


text:
Can MSstatsTMT be used for both MS2 and MS3 runs?
MSstatsTMT is an R library for statistical analysis of TMT labeled mass-spectrometry (MS) proteomics data. The documentation does not explicitly mention the kind of MS that can be processed. The examples contain comments that contain both MS2 and MS3, but I would like to make sure that is indeed the case. Can anyone confirm that and maybe provide an example script?


text:
Discordance in gene signature behavior between bulk and single-cell RNASeq
The objective of the following analysis is to identify an activation signature of a specific phenotype on bulk RNASeq and to apply it to single-cell RNA-Seq, in order to identify the population of activated cells in scRNA-Seq samples.
In the first step I used DESeq2 to identify genes that increase during different timepoints (I have 4 timepoints):

Bulk samples are normalized using VST (Variance Stabilizing Transformation - DESeq2)
For each timepoint, differentially expressed genes are computed with respect to the previous timepoint (i.e. 2 vs 1, 3 vs 2, 4 vs 3)
Differentially expressed genes that increase with timepoint (log-FC > 1, padj < .05) are selected (n=24)

As a control, I plotted the activation score (computed with scanpy) based on the 24-gene signature applied to the bulk samples from which it has been extracted. As you can see, the score always increases from the first timepoint (blue) to the last (red).

However, when I use the same scoring procedure on a single-cell RNASeq sample containing the same populations, I think something may be wrong.
In the following the violin plot of the extracted activation score on an unsupervised clustering on single-cell data. Excluding cluster 4, which we identified as contamination, while cluster 1 and cluster 0 are supposed to have the correct degree of activation, clusters 2 and 3 have a lower degree of activation with respect to cluster 0 (which we identified as the first timepoint - lowest activation on bulk).

What may be the causes of this behavior? Right now the possible explanation I came out with is that:

Clusters 2 and 3 have global shifts on gene expression that make TPM normalization over-normalize activation genes

Are there other possible explanations for this behavior, and how to check for them?

How do you know the temporal relationships between clusters in the single-cell data?

@TimStuart, I correlated the expression profile of each cell in each cluster with the expression profile from independent single-cell samples in which I know the temporal relationship. Cluster 0 is more correlated to the first timepoint, while cluster 1-3 are more correlated to the consecutive timepoints.

You assume that your gene signature is originating from a single cell type, from the "activated cells". This so called activation can result from an interplay between different cell types (this is just a speculation as I do not know your experimental settings and the biology here). A heatmap of genes vs clusters would show the "source" of the signal. In my field (cancer), we have observed that tightly related coexpressed gene modules (bulk RNA-seq) might originate from diverse cell types (scRNA-seq).

One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).
Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).
The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.


text:
Mix globbing and wildcards when specifying rule input
Consider the following scenario.
├── sample-alice
│   ├── sequence_1.fastq
│   ├── sequence_2.fastq
│   ├── ...
│   └── sequence_n.fastq
└── sample-bob
    ├── sequence_1.fastq
    ├── sequence_2.fastq
    ├── ...
    └── sequence_m.fastq

I'm trying to write a single rule in my Snakefile that will preprocess each sample generically. The rule will execute a single command that takes all of Fastq files associated with the sample simultaneously. The samples do not have the same number of Fastq files.
If it was a single sample, I could do something like this.
rule preprocess:
    input: glob('sample-alice/*.fastq')
    output: 'alice-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

On the other hand, if there was a single Fastq file for each sample, I could use a wildcard to write a generic rule like this.
rule preprocess:
    input: 'sample-{samp}/sequence_1.fastq'
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

Is it possible to combine these two approaches to write a rule for handling multi-file samples generically?

I'm not sure if this exactly solves your question, but in one of my workflows I use the directory as input (as specified in the config.yaml). The shell command then uses {input}/*.fastq.gz See some examples here: https://github.com/wdecoster/nano-snakemake/blob/master/rules/align.smk

A nice solution comes from a colleague of mine, as seen in this Snakemake workflow. The trick is to access the wildcards programmatically using an anonymous (lambda) function. In my example above, it would be implemented as follows.
from glob import glob

rule preprocess:
    input: lambda wildcards: glob('sample-{samp}/*.fastq'.format(samp=wildcards.samp))
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

Note that your Snakefile needs a from glob import glob for this to work

And it's a bit more succinct using fstrings: glob(f'sample-{wildcards.samp}/*.fastq')


text:
Differential expression analysis when nested effects
I have 3 tumour samples from 3 patients from an experiment, I also downloaded 10 normal samples from TCGA. My design is like this
> mycols
               condition batch
normal1           normal     1
normal2           normal     1
normal3           normal     1
normal4           normal     1
normal5           normal     1
normal6           normal     1
normal7           normal     1
normal8           normal     1
normal9           normal     1
normal10          normal     1
OESO_036_a_RNA    cancer     2
OESO_013_a_RNA    cancer     2
OESO_005_a_RNA    cancer     2
>

But DESeq2 returns error when trying for differential expression analysis
> dds <- DESeqDataSetFromMatrix(countData=counts, colData=mycols, design=~batch+condition)
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')
In addition: Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

Then how I could extract deferentially expressed genes between tumour and normal samples while they are coming from two different batches?
Any help please?

Did you read the vignette section "Model matrix not full rank"? The problem here is that there is no way you can distinguish the effect of being normal and the effect of being sequenced on batch 1, you have the two effects affecting the same samples!

There should be a solution for this case , is not it?

Not unless you have some normal and cancer samples that you can resequence in a third batch. (But perhaps you don't have a batch effect, but I wouldn't bet anything on this)

llrs is right.  "Matrix Model not full rank" means that two of your conditions are confounded.  You can't get your normals from one place and your experimentals from another
And no, there generally is no clever way out of this.  For all you know, any differences you see between your normals and your cancers are 100% due to being sequenced somewhere else, and have nothing to do with real biology.
If you want to run this through DESeq, drop batch effect.  But tell whoever you are reporting the results to that batch is hopelessly confounded with cancer state, and you can't tell which differences are due to what.

Sorry I was designing an expensive RNA-seq experiment so I must be wise before invest any money; could I please ask you if the design gives what I want? You please imagine I have 1- fibroblast surrounding the tumor, 2- tumor itself and 3 - I have co-cultured  fibroblast  with tumor; My ultimate goal is getting genes coming from the interaction of fibroblast and tumor, so I am going to use DESeq2 to for differential expression the condition would be co-culture vs fibroblast + tumor; Do you think this is right?


text:
What is the recognized way to reference accession data (.sra files) from NCBI as a URI?
I would like to assign a URL/URI to accession runs (i.e. SRR IDs) from the NCBI short read archive.
In program listings, I used to list them as ftp:// URLs, because they were available on their servers as such, e.g. SRR3418242:
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR341/SRR3418242/SRR3418242.sra

But this approach is deprecated, plus, not all runs are available this way anymore. The now recommended approach is to use sra-tools to download the data. This tool takes the SRR Ids directly, i.e. prefetch "SRR3418242", and resolves them to one of many possible locations (home directory, http, etc).
Is there an official URI scheme that can be used to reference accession runs from the short-read-archive ? I would like to disambiguate them from other types of files and URLs (S3 keys, http, file). For file manifests, spreadsheets, and publications, is there something similar to srr://SRR3418242?

No, there isn't something like srr:// and there won't be. ENA via aspera is almost always the preferred way to download data as long as ENA has it.

Thanks. My concern is the namespace. I'd like to know a nomenclature which doesn't prescribe a particular software or particular read archive. As long as run numbers and sequence data agree on both archives...

The nomenclature is no longer being followed as NCBI is moving all the files to S3/Google cloud. pysradb allows you to fetch the metadata and get the URLs for each individual SRR. See examples in this notebook.
You can also download an entire project by doing:
pysradb download -p <SRP>


text:
"Bad substitution" error while trimming adapters
I have this script to trim adaptors from my fastq sequencing files
for PAIR in $(fastq | sed 's/_R.//g' | uniq)

I am gitting this error: bad substitution
Could you help

Welcome to Bioinformatics.StackExchange! I'd guess this error message is being generated by the sed command. But it looks like there are several potential problems with your shell command, the braces around variable names and the listing of the input files (the fastq | seq | uniq sequence) in particular. It doesn't look quite right, but to help you it would be helpful if you could give a bit more details about the names and locations of your input files. If there are many, 1 or 2 examples would suffice.

Hello again, Momo. Your edits actually made it harder for us to answer the question rather than easier. I've voted to close this question, but I'll revoke that vote if you can provide a bit more detail as discussed in my first comment.


text:
How to extract gene expression tables from this GEO dataset?
I've downloaded this GSE43013 dataset using GEOquery in R. My understanding is that it contains expression data from liver, kidney, and brain for multiple species. I would like to produce gene expression tables for each species and tissue type. The goal would then be to perform PCA and regression analyses to see if there is a relationship between expression of certain genes and e.g. taxonomic position, diet, etc.
However, I have never worked with gene expression data before so I have very little understanding of the GSE file. I have failed to extract any kind of expression data from the file. Solutions that I have found online do not seem to be working for this dataset (for example producing empty tables where I assume there would normally be data on expression levels). Because of my inexperience I cannot determine why they are not working in this situation.
I would appreciate some help in understanding the format of this file, but most importantly, in figuring out how to pull expression data from it, if that is possible.
EDIT: I still haven't managed to pull data from the GSE file, but looking through the samples associated with the data set (such as this one), I can see that there are files formatted like this:
mus.lv.3
ENSMUSG00000000001  26.3376
ENSMUSG00000000003  0
ENSMUSG00000000028  0.501533
ENSMUSG00000000031  0.116078
I assume that the values in the left column refer to specific genes. Perhaps the values in the right column represent expression intensity? If so, would these values be normalized and comparable between species as they are?

Possible duplicate of Extracting expression data from GSE dataset downloaded from GEO

Assuming that you have downloaded processed file, according to your link, it is either plain text (compressed) or GTF. If you would have plain text you probably would have posted a snapshot or sample of it, so I assume you have the GTF. According to format description for GTF-files:

Format consists of one line per feature, each containing 9 columns of data, plus optional track definition lines.
The 9 columns (or fields) are tab-separated and all fields contain a value (missing/NA-values are usually denoted with ".").
Fields are:

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

You can find more information on it here and here.
IF you don´t have the file and just have the query output, you probably will find the answer here, where extraction of expression data has been discussed in more detail.


text:
Kinetics "Master Equation"
In stochastic (probabilistic) models, the evolution of a system over time is calculated using a "Master Equation" which also factors the external and internal noise / perturbations in order to more accurately replicate the behavior of a biological system (e.g. cell growth). I am unable to understand what is it, specifically, about the "Master Equation" that allows computational biologists to replicate a biological network with noise?
The Master Equation is the phenominological description of a non-stationary (i.e., evolving with time) process.

Welcome to the site. Perhaps it would be easier to answer if you pointed towards a paper or a course book where you found this reference. I have never heard of this Master Equation (but maybe this is specific to a sub-field or something)

@mobs By "noise" here are you referring to the probabilistic change in state of individual members? Regarding perturbations, if you mean external perturbations then they're not always accounted for (they would need to be made explicit in the differential equations).


text:
RNA-seq analysis of mixed viral/host reads with salmon
I did RNA-seq of mammalian cells infected with pox virus. Now, I have read files which contain both host and virus reads. I want to align the reads both to host and viral genome. I was thinking I could concatenate the host and virus genome into one file and run Salmon against this concatenated genome.
However, salmon recommends transcriptome file for assembly which are not available for viruses. Virus genome are available as GenBank or GFF3 format from NCBI. Is there any way I can concatenate these formats into the format that can be used by salmon? Or is there any way around to use virus genome as reference in salmon?

Citing the website:
Salmon is a tool for quantifying the expression of transcripts using RNA-seq data. Salmon uses new algorithms (specifically, coupling the concept of quasi-mapping with a two-phase inference procedure) to provide accurate expression estimates very quickly (i.e. wicked-fast) and while using little memory.
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods.

To be clear, salmon does not align reads to the genome or transcriptome per se, but uses k-mers to determine which transcript from which a read is likely derived. If you really need the alignments, you'll need to use a different tool. But if you're just trying to do transcript quantification, tools like salmon or kallisto are MUCH faster and just as accurate as alignment-based tools.

You provide Salmon with a transcriptome fasta... so merging the human transcriptome with the pox virus genome fasta file should work. You don't mention that in your post but the viral genome should exist as a fasta. You can extract the sequences from genbank or maybe you can find our virus here: https://www.ebi.ac.uk/genomes/virus.html
You can find the coding transcriptome at the link below in the cdna folder.
The non-coding transcriptome is in the ncrna. You should probably merge them:
ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/

Yes, I have virus genome in fasta format. I was wondering if I just concatenate in plain fasta format , the features will be missed. How do I merge so that I will have virus features included in the merged file ?  Thank you.

If you want the virus features then you would need to find the ORF/CDS regions and subset your original fasta. You could try something like: https://www.ncbi.nlm.nih.gov/orffinder/ if they are not annotated. Alternatively they might already be present in Genbank/NCBI so you could just download them from there - although I'm not familiar enough with this to tell you how to do it in a high-throughput manner.


text:
Batch alignment of inconsistently named Fastq files
I have numerous gzip-compressed paired-end Fastq files with ChIP-seq data that I would like to align with bowtie2. I confirmed that bowtie2 takes .gz files, although I do see this warning when I run for two samples:

Warning: gzbuffer added in zlib v1.2.3.5. Unable to change buffer size
  from default of 8192.

My files are named like so.
10_S10_L002_R1_001.fastq.gz  12_S12_L002_R1_001.fastq.gz  2_S2_L001_R1_001.fastq.gz  4_S4_L001_R1_001.fastq-003.gz  6_S6_L001_R1_001.fastq-002.gz  8_S8_L002_R1_001.fastq.gz
10_S10_L002_R2_001.fastq.gz  12_S12_L002_R2_001.fastq.gz  2_S2_L001_R2_001.fastq.gz  4_S4_L001_R2_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz  8_S8_L002_R2_001.fastq.gz
11_S11_L002_R1_001.fastq.gz  1_S1_L001_R1_001.fastq.gz    3_S3_L001_R1_001.fastq.gz  5_S5_L001_R1_001.fastq-002.gz  7_S7_L002_R1_001.fastq.gz      9_S9_L002_R1_001.fastq.gz
11_S11_L002_R2_001.fastq.gz  1_S1_L001_R2_001.fastq.gz    3_S3_L001_R2_001.fastq.gz  5_S5_L001_R2_001.fastq-003.gz  7_S7_L002_R2_001.fastq.gz      9_S9_L002_R2_001.fastq.gz

The bash commands I tried are below. When I run this test command:
for i in $(ls *.fastq*.gz | grep -v 6 | rev | cut -c 16- | rev | uniq); do
    echo $i;
done

...I get this output.
12_S12_L002_R1_001.fastq.gz  12_S12_L002_R2_001.fastq.gz

When I run this:
for i in $(ls *.fastq*.gz | rev | cut -c 16- | rev | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}R1_001.fastq.gz -2 ${i}R2_001.fastq.gz | samtools view -bS - >${i%}.bam
done

I do see one BAM for my paired end files. But as I see the second set of files which are named:
6_S6_L001_R1_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz

if I run the same script over them, I would get an error as the loop can't read the paired end file.
So the tedious way is I can make a separate folder for these sets of data with fastq-002.gz and fastq-003.gz and can run the script by tweaking the cut parameter.
But I would like to know what is the programmatic way of resolving the issue.

I strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud.

Can't you just rename the files?

@terdon that's how our collaborator gave us and i dont want to mess up with the file name yes i can if i dont find any other way..

" strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud." i read a bit about it never tried it , i need to learn and give it a try

It's good practice to treat files from collaborators as read-only, so not changing the filenames is a good idea. But if the files are inconsistently named, it's a very good idea to create consistently named symlinks to those files. Well organized filenames and directory structures will make it SO NICE AND EASY to automate batch data processing.

yes i did that finally

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to:
for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz |
        samtools view -bS - >${i%}.bam"
done


text:
What is the difference between a MeSH concept and a MeSH descriptor
Medical Subject Headings (MeSH) is a very comprehensive vocabulary to index the biomedical literature on Medline or PubMed (which is the major Biomedical literature database). MeSH is arranged into Descriptor > Concept > term hierarchy, but I still do not understand the difference between a MeSH descriptor and a MeSH concept.
For example, for Knee

MeSH descriptor is D007717
MeSH concept is M0012066 (Under the concepts section in the link)

Even after reading the description here and here, I am still not sure about the purpose of MeSH concept and MeSH descriptor.
Could anyone please shed some light on it?

MeSH is organized into Descriptor > Concept > Term hierarchy, where all the synonymous <term>'s are grouped together under a <concept>. For example, consider the MeSH terms,

Blood-Brain Barrier
Hemato-Encephalic Barrier
Brain-Blood Barrier

All these terms are group under the umbrella concept of "Blood-brain barrier".


text:
Storing FindAllMarkers results in Seurat object
I am currently working on multiple datasets where each is managed by a separate Seurat object. For each I'd like to also compute the marker genes using the FindAllMarkers function. Where would it be a good practice to save the markers returned from the FindAllMarkers function? I tried adding an extra slot to the Seurat object with no success.
I saw the calc.params slot and though I could add an item to that list. However, it doesn't seem like the right place as that slot is a list of the parameters used for the different functions executed on the data.
Currently I have some wrapping object that has a slot for the Seurat object and a slot for the results from the FindAllMarkers. Is there any better place to put it? Somewhere within the Seurat object?

Welcome to the site Gilad. In general it is discouraged to directly access and modify the slots of an object directly. If the object itself doesn't add the information it might be a good reason to not attempt to do it

@llrs - True which is why I'm checking to see if there is any suitable slot that is used for this kind of extensibility. As in the example of calc.params where it is designed in such a way where it is ok to add items to that list

You can stash anything you like in the misc slot (present in both v2 and v3 Seurat objects). I often use it for storing marker information to help with organization.


text:
Downloading Fasta for each chromosome of a complete genome
Whatever I am googling I am not able to find fasta file of each chromosome for hs37d5.fa. hs37d5.fa whole fasta file however is here
ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
Do you know where I can get fasta file for each chromosome for this genome?

Is there something stopping you from downloading the entire FASTA and then splitting it? Also, there are quite a few typos in your post making it a little difficult to understand - can you fix those please? Thank you!

@oka I rolled back your edit because it seemed to be introducing almost as many gramatical issues as it was fixing.

@FereshTeh Sorry, but does your comment have anything to do with the question at hand?

Sorry you right, I deleted my comment

Personally, I would download the entire file and split it using something like fastaexplode from the exonerate package. On an Ubuntu machine, you can do that with:
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
gunzip hs37d5.fa.gz
sudo apt install exonerate
fastaexplode hs37d5.fa

Alternatively, you could do this on the fly while downloading:
wget -O - ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz |
  gunzip | awk '/>/{name=$1; sub(/>/,"",name)} {print >> name".fa"}'

Thank you  very much, [fi1d18@cyan01 ~]$ sudo apt install exonerate
[sudo] password for fi1d18:     So I am not able to install from root, your alternative suggestion gave me all chr fasta, but why you are saying on fly?

@FereshTeh yes, you need root access to install software. If you are not running this on your machine, then next time remember to mention that in your question. As for "on the fly", that's because that command downloads the file to standard output, decompresses it as it is coming in and then writes the output to each  .fa file, so it is creating the files on the fly.

It is common that reference genomes will be published only as a complete genome, without separate files for each chromosome. Fortunately, a google search for "split fasta by chromosome" will give you many options for downloading a complete genome Fasta file and splitting it by chromosome. The top answer on this biostar thread looks promising.
pip install pyfaidx
pyfaidx -x sequences.fa

Alternatively, I typically use a homebrew Python script. If there are unlocalized/unplaced scaffolds in the file, this will create a bunch of files with ugly unexpected names. These are easy to ignore or delete.
# usage: ./split.py < genome.fasta
from Bio import SeqIO
import sys
reader = SeqIO.parse(sys.stdin, 'fasta')
for record in reader:
    outfile = 'chr{chrid}.fasta'.format(chrid=record.id)
    with open(outfile, 'w') as outstream:
        count = SeqIO.write(record, outstream, 'fasta')

With a little bit of searching, you will find a wealth of helpful information on blogs, forums, and library documentation.