First off,

Don’t use RPKMs.

They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use FPKMs, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count fragments).

Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:

$$
\text{TPM}_\color{orchid}i =
    {\color{dodgerblue}{\frac{x_\color{orchid}i}{{l_\text{eff}}_\color{orchid}i}}}
    \cdot
    \frac{1}{\sum_\color{tomato}j \color{dodgerblue}{\frac{x_\color{tomato}j}{{l_\text{eff}}_\color{tomato}j}}}
    \cdot
    \color{darkcyan}{10^6}
$$

where

$\color{orchid}i$: transcript index,
$x_i$: transcript raw count,
$\color{tomato}j$ iterates over all (known) transcripts,
$\color{dodgerblue}{\frac{x_k}{{l_\text{eff}}_k}}$: rate of fragment coverage per nucleobase ($l_\text{eff}$ being the effective length),
$\color{darkcyan}{10^6}$: scaling factor (= “per millions”).

That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:

fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}

Here, the effective length is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.

This function handles one library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:

tidy_expression = tidy_expression %>%
    group_by(Sample) %>%
    mutate(FPKM = fpkm(Count, col_data$Lengths))

However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:

tidy_expression = expression %>%
    gather(Sample, Count)

This equation fails if all your counts are zero; instead of zeros you will get a vector of NaNs. You might want to account for that.

And I mentioned that TPMs are superior, so here’s their function as well:

tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}


How do I infer someone's haplotype of a gene from data on several of the gene's SNPs? For a a few SNPs from given gene, I know there happen to be several haplotypes that most of the population has. But from a person's data on those SNPs, we have two alleles at each SNP, so it's not trivial to match the SNP data to infer the person's haplotypes. When I search this I find tools that seem to be for haplotype reconstruction and SNP imputation, such as SHAPEIT, but this is not what I'm looking for.
Might someone know of a tool for this inference of known haplotypes from SNP data? It would be helpful even if someone could give me a name for this problem, as I'm not sure of the right term to search.

It depends on the scale of the data (i.e. number of individuals, number of SNPs). Generally this can be done by matching up individual SNP genotypes to haplotype tag SNPs. The specific combination of genotyped SNPs is used to determine the underlying haplotype, based on tag SNP frequencies for each haplotype.
More details of the process of identifying tag SNPs can be found here:
https://www.broadinstitute.org/haploview/tagger


Protein sequence from patient data Currently, I am working on NGS data and my aim is to get significance prediction of variants present in the vcf file. As we know about SIFT Score for significance score prediction, I am trying to understand how this score works. When I read its reference paper, in the method section, the procedure to determine the sift score starts with patient protein query sequence.

My question is how to get patient protein sequence from patient data (either vcf file or bam file) (which is not explained in the paper). Rest I can follow by other algorithms like PSI-BLAST & MOTIF.

Any suggestions....

Thanks

Ensembl now have a tool called Haplosaurus (still in beta) which will convert your phased VCF into actual protein sequences.

  Haplosaurus

  haplo is a local tool implementation of the same functionality that
  powers the Ensembl transcript haplotypes view. It takes phased
  genotypes from a VCF and constructs a pair of haplotype sequences for
  each overlapped transcript; these sequences are also translated into
  predicted protein haplotype sequences. Each variant haplotype sequence
  is aligned and compared to the reference, and an HGVS-like name is
  constructed representing its differences to the reference.

  This approach offers an advantage over VEP's analysis, which treats
  each input variant independently. By considering the combined change
  contributed by all the variant alleles across a transcript, the
  compound effects the variants may have are correctly accounted for.

  haplo shares much of the same command line functionality with vep, and
  can use VEP caches, Ensembl databases, GFF and GTF files as sources of
  transcript data; all vep command line flags relating to this
  functionality work the same with haplo.

  Usage

  Input data must be a VCF containing phased genotype data for at
  least one individual; no other formats are currently supported.

  When using a VEP cache as the source of transcript annotation, the
  first time you run haplo with a particular cache it will spend some
  time scanning transcript locations in the cache.

./haplo -i input.vcf -o out.txt -cache


How to access the patient-derived xenografts (PDXs) repository? I saw this nature news item, it sounds that Cellminer is obsolete, is it right?

What are the new tools to analyze the "new cell lines"? Where's the PDX repository?

What about the European initiative that the article refers to?

I would be happy to see correlation between drugs and genes and similar capabilities as we can see on Cellminer.

I already see for links when I try to search. Take a look at them below. Also if you can be precise about kind of data are you trying to find, it would be better

https://pdmr.cancer.gov/

http://cdt.northwestern.edu/news/patient-derived-xenograft-repository-now-available-researchers

http://www.epo-berlin.com/epo-tumor-models-xenografts.html

http://www.proxe.org/

http://data-analysis.charite.de/care/

http://www.oasis-genomics.org/

http://www.cbioportal.org/

Hope one of these above will serve the purpose for you.

Update: added some more links that might serve. As for the Cellminer the last update was 2012.


Demultiplexing FASTQ file without index information I am trying to understand how data that was uploaded to SRA (https://www.ncbi.nlm.nih.gov/sra?LinkName=biosample_sra&from_uid=4510743) can be analyzed with the assumption that the FASTQ file should contain replicate information (3 replicates within the sample). I downloaded the file and ran it with Nextflow, and it was mapped to the human genome and determined to contain unstranded reads of 50bp, both directions.
Which steps can I take to infer if the FASTQ file indeed is not demultiplexed?
Many thanks,
Geo
The header of the FASTQ looks like this:
@SRR3467208.1 1/1
NTCTCTATGTCCACTCTGGAGCCTTTAAGTGCCACTTGAGGGCCTTTAAC
+
#4:D;DDEHHFHHIJJIJJGHGHJJI@HFIIJJHJJIG?9CFHIIIHHHG
@SRR3467208.2 2/1
NGTGCTGTGTGTGCATGTGTGTGCGTGTGTGTGCTGTGCGTTTGTGTGTG
+
#1:B:BDDFFFFFFFBGFHIFIIFGFFGFGFGIIIIIIEFGGIIFFIFII
@SRR3467208.3 3/1
NGTCACTTGGCATCTGGCCATCGGGCTGGATGCCGTGTTCCAGGCAGTAG
+
#1:A=DDDB<<CFAFGI@FGG9A@68CEG;DB?<DGGGE49?BBF2B=FF
@SRR3467208.4 4/1
NTCCGGGGCTCCAGCAACCAGAAGAGGGAAAAGCTGTCTTCGGTGATGTC
+

I would suggest looking at sample metadata:

Library:

Name: TruSeq RNA

Instrument: Illumina HiSeq 2000

Strategy: RNA-Seq

Source: TRANSCRIPTOMIC

Selection: PolyA

Layout: PAIRED

You can see that it's paired. You can look at the spot descriptor at the page you link to see more information about the design.
There is no multiplexing information available there. I would instead look at the overall project and its metadata to find the different replicates. It looks like there are 2 samples x 3 replicates each, deposited independently. So it is likely to already be demultiplexed, and you can just download them individually.
"Unstranded" isn't really a meaningful term for this, if I understand correctly. Forward and reverse read pairs don't correspond to forward and reverse on the chromosome, but to the way that sequencing-by-synthesis happens on the machine. So the mapping orientation is not really relevant.
I infer that you are working from a single FASTQ file. If you are working from a single FASTQ file for paired-end out of SRA, then the FASTQ file is not yet resolved into F+R reads. I would suggest rerunning the FASTQ conversion with --split-3 or similar option.


Will blastn remove sequences from a search with low identity? I'm using this command (excuse the duplicate naming, I know it's bad form):

blastn -query cm-seqs/combined_seqs.fna -db combined_seqs.fna -out cm-matched.txt -num_alignments 1 -outfmt 10

to take a set of fasta sequences to blast against another data base I built with blast.

My query file has 4,364,417 sequences. The files resulting from the above command spat out 4,362,639 results. The difference is about 1800. Is this due to poor alignments for some sequences that didn't pass any of blasts default parameters? I assumed the result would be 1 "best" match in the database for each query sequence.

Yes, if a query sequence has no good alignment, blast will not return a result for it. Your command will indeed output only the best hit for each query sequence, but there is no reason to assume all of your sequences will have a match.

Those that don't have a match won't be included in the output. This sounds perfectly normal unless you are blasting the same input file against a DB created from that file. Otherwise, you don't expect all of your sequences to necessarily have a match.

Also note that since the default blastn settings are quite permissive (usually an e-value of 10, for example), don't make the mistake of assuming that a blast hit means anything. The details will always depend on the specific search you are running (I have discussed this a little in my answer here), but in the vast majority of cases, any hit with an e-value > 1 can safely be ignored. In most cases, you are only interested in e-values well under 1, actually. So you probably already have all sorts of spurious results in there anyway.

For more details, please edit your question and explain what the query and DB files are.


How to convert GFF3 to GTF2 I would like to convert a file in gff3 format to a gtf2.2 format.

The reason why I would like to do this is: I have a set of transcripts assembled by a bunch of different software (and using RNA-seq data from different sequencing technologies) and I would like to compare each of those sets with the reference transcriptome annotation for that species (D. melanogaster). I already asked for the community advice about how to proceed on that, but to run the suggested software (TACO and cuffmerge) I need to have a GTF file containing the experimental transcripts to compare to the reference transcripts (GTF).

Up to now, I unsuccessfully tried:

gffread utility in the Cufflinks package (gffread input.gff3 -T -o output.gtf): this results in an empty output.gtf file and an empty log file (used Cufflinks v.2.2.1) - I contacted the authors via their Google group but haven't heard of them yet
gff3_to_gtf utility in the GenomeTools (gt) package (gt gff3_to_gtf input.gff3 -o output.gtf): the output is not created and the log file is not informative - I contacted the authors via their mailing list but haven't heard of them yet
GFF3_to_GTF utility in the FML package (./gff3_to_gtf_converter.pl input.gff3 output.gtf): the output just contains a header (##gff-version 2.5) and the log is empty

The gff3 file was created as output of GMAP, and contains the transcripts as found by alignment to the reference (specifying option -f gff3_match_cdna).

[Edit: what I found out a poteriori, is that this format is not a standard gff3, thus the conversion is not trivial...]

Here is the beginning of the gff3 file I tried to convert:

$ head r9_gmap_match-cdna.gff
##gff-version   3
# Generated by GMAP version 2017-04-24 using call:  /home/aechchik/software/gmap-2017-04-24/bin/gmap.sse42 -d gmapidx -D /scratch/beegfs/monthly/aechchik/isoforms/ref/chromosomes/ -f gff3_match_cdna -n 0 -t 20 r9_2d.fasta
2L    gmapidx    cDNA_match    18442664    18443024    79    -    . ID=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1;Name=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d;Target=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d 141 489;Gap=M13 I10 M8 D3 M8 D3 M32 D1 M4 D2 M5 D3 M34 I3 M6 D1 M7 D1 M5 D1 M30 D2 M16 I1 M8 D3 M10 D1 M5 D1 M6 I1 M8 I1 M8 D1 M33 D3 M33 I1 M28 D3 M25
###
3R    gmapidx    cDNA_match    15853880    15855465    96    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 80 1645;Gap=M46 D2 M12 D2 M157 D3 M4 D1 M50 I2 M3 I1 M12 I1 M68 I1 M66 I1 M53 D2 M47 D1 M16 D1 M35 D1 M155 D1 M28 D1 M166 D2 M47 D1 M8 D4 M69 D1 M28 D1 M5 D1 M12 D1 M202 I1 M115 D1 M61 I2 M7 D1 M7 I1 M36 D2 M41
3R    gmapidx    cDNA_match    15855529    15855742    97    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 1646 1856;Gap=M21 D1 M68 D1 M85 D1 M37
###
X    gmapidx    cDNA_match    14837810    14838142    93    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 74 406;Gap=M13 D1 M182 I1 M14 I2 M56 I1 M30 D2 M21 D1 M13
X    gmapidx    cDNA_match    14837470    14837753    92    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 407 688;Gap=I2 M5 I1 M64 I1 M44 D1 M9 D5 M31 D3 M23 I1 M19 I1 M25 I1 M26 D1 M19 I1 M9
###

The following bit of python code should work:

#!/usr/bin/env python
import sys

lastTranscript = [None, None, None, []]  # ID, chrom, strand, [(start, end, score), ...]

def getID(s):
    """Parse out the ID attribute"""
    s = s.split(";")
    for k in s:
        if k.startswith("ID="):
            return k[3:]
    return None

def dumpLastTranscript():
    """Print the last transcript"""
    bounds = sorted(lastTranscript[3])
    print("{}\tgmapidx\tgene\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    print("{}\tgmapidx\ttranscript\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    for start, end, score in bounds:
        print("{}\tgmapidx\texon\t{}\t{}\t{}\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], start, end, score, lastTranscript[2], lastTranscript[0], lastTranscript[0]))

def handleLine(cols):
    """Handle a single line, appending the exon bounds to the previous if relevant"""
    ID = getID(cols[8])
    assert(ID is not None)
    if lastTranscript[0] is not None and lastTranscript[0] != ID:
        dumpLastTranscript()
        lastTranscript[3] = []
    lastTranscript[0] = ID
    lastTranscript[1] = cols[0]
    lastTranscript[2] = cols[6]
    lastTranscript[3].append((int(cols[3]), int(cols[4]), cols[5]))

f = open(sys.argv[1])
for line in f:
    if line.startswith("#"):
        continue
    cols = line.strip().split("\t")
    handleLine(cols)

dumpLastTranscript()
f.close()

If you saved that as gff2gtf.py and made it executable, the usage would be ./gff2gtf.py foo.gff3 > foo.gtf. With the example that you provided in your post, the result is:

2L  gmapidx gene    18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx transcript  18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx exon    18442664    18443024    79  -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
3R  gmapidx gene    15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx transcript  15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15853880    15855465    96  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15855529    15855742    97  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
X   gmapidx gene    14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx transcript  14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837470    14837753    92  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837810    14838142    93  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";

That at least looks like correct gtf 2.2 format to me.


confidence ellipses for MDS plot in edgeR? Is it possible to draw e.g. 95% confidence ellipses around samples from the same group on the results from the plotMDS function under edgeR? If so, how?

I can see that there would be a three step process to doing this:

Merge counts from all samples in the group and then resample pseudo-replicates from this. If x is a matrix with samples in the group being columns and genes being rows

s <- sample(row.names(x), n = mean(colSums(x)), probs=rowSums(x)/sum(rowSums(x))
stab <- table(s)
s <- as.vector(stab)
names(s) <- names(stab)

Do this thousands of times.
You would then want to project those onto your MDS plot space - i'm less clear about how you could do this without perturbing the space itself.
Calculate the parameter for an ellipise that would contain 95% of these points. Again, I'm not so sure about how to do this, but I think the car pacakge might be a good place to start looking.


Specific cell type identification in Single Cell Sequencing In order to define which cell is of which type we need to identify a set of rules, for instance neurons should express one of the following: Thy1, Rbfox3, MAP2, Camk2b, Gad1,Cck, Reln, and should not express any of the following: cd45, Tmem119, CD11b, ... and others. Maybe Rbfox3 should always be expressed in any neurons. So, these rules need to be compiled manually as I understand to identify each cell type, from the literature. I would assume that many people already faced this issue, and maybe came up with software to do that? Are there any software where we would be able to supply a matrix of single cell expression data and have the types assigned with certain degree of certainty?

I know that maybe we can also use heatmaps for that, like:

cell type specific genes - heatmap using rank based approach

But again, this seems too manual work. Besides, I do not think that we can assign different neuronal types based on just a heatmap.

I do not know of such software.

However, I believe this effort is a bit misdirected. The purpose of single-cell sequencing is to get a better understanding of cells; their heterogeneity and functional diversity or developmental / biological processes such as differentiation, using a higher "resolution" method. In other words, if we had the methods you are asking for, there would be no need for the single-cell experiment.

Running your data through some pipeline utilizing previous knowledge has the danger of forcing ideas on the data, rather than seeing what the data tells you: it would be better to try and understand and explain the biologically relevant heterogeneity and diversity in your cells, together with a critical comparison of cell type characteristics (expression of genes, pathways) with knowledge described in literature.
This question is also very similar to the one you previously asked and was answered.

Also good to recognize that your approach has a lot of implicit assumptions, such as: cells can be clearly categorized as +/- expression for each gene, or it ignores the systems level (networks and pathways), thus again forcing concepts on the data.

This said, one method would be to construct a table listing cell type markers as described in the linked answer above, then write a script that determines a cutoff value for expression of genes in your data (see this), then ranks cell types for each cell. For example, you can measure the number or proportion of uniquely genes expressed (although that's very simplistic).

Another option, if you are familiar with machine learning, is to train a classifier on an annotated dataset, and then use that on new data.

Also see a convenience function below, which requires a table of marker genes for each cell cluster (i.e. (the Seurat::FindAllMarkers() output), and a reference df listing genes and corresponding cell types in its HGNC_symbol and Cell_type columns, and returns the table with listing clusters, their marker genes and corresponding cell types.

getCelltypes <- function(markers, reference) {
  marker.celltype <- markers
  marker.celltype$Cell_type <- marker.celltype$gene
  marker.celltype$Cell_type <- with(reference, celltype[match(marker.celltype$Cell_type, gene)])
  return(marker.celltype)
}

Many variants of this reference table and function can be created, feel free to modify it.

In the future, I think we can expect such classifier software that uses reference data from the Human Protein Atlas, HCA and similar projects.


What is the minimum number of individuals per group required for ADMIXTURE analysis? What is the minimum number of individuals required per group (population/species) for an ADMIXTURE analysis to be appropriate?

Is there a way to determine this number mathematically based on the amount of genetic data I have or other criteria?

Are there any citations I can use to figure this information out?

The following paper provides an in-depth exploration of factors affecting power of admixture analyses under a number of different models. Although the authors state they "provide the open-source R software AdmixPower, to perform power and sample size analysis", in fact it seems the software is not available.

AdmixPower: Statistical Power and Sample Size Estimation for Mapping Genetic Loci in Admixed Populations

Tool availability:

https://research.cchmc.org/mershalab/Tools.html

Site asks for a password, apparently there is no option to sign up.


Pattern matching and removing from the name colnames(countdata) = gsub("Aligned.sortedByCoord.out.bam_rem.bam", "", colnames(countdata))
colnames(countdata)
dput(colnames(countdata))
    c("BAF60A_KD_1_S13", "BAF60A_KD_2_S14", "BAF60A_RA_1_S17", "BAF60A_RA_2_S18",
    "BAF60A_VD3_1_S15", "BAF60A_VD3_2_S16", "BAF60C_KD_1_UI_S1",
    "BAF60C_KD_2_UI_S2", "BAF60C_KD_RA_1_S11", "BAF60C_KD_RA_2_S12",
    "BAF60C_KD_VD3_1_S9", "BAF60C_KD_VD3_2_S10", "VC_1_UI_S3", "VC_2_UI_S4",
    "VC_RA1_S5", "VC_RA_2_S6", "VC_VD3_1_S7", "VC_VD3_2_S8")

After removing the above pattern I want to remove this part from each of the colnames _S[anynumber] so that I can end up like this as an example  BAF60A_KD_1 but not sure how to make that pattern and match the same and remove or replace it with empty space. I tried as below which ends up like all the digits are removed.

colnames(countdata) = gsub("[[:digit:]_S]", "", colnames(countdata))

colnames(countdata)

c("BAFAKD", "BAFAKD", "BAFARA", "BAFARA", "BAFAVD", "BAFAVD",
"BAFCKDUI", "BAFCKDUI", "BAFCKDRA", "BAFCKDRA", "BAFCKDVD", "BAFCKDVD",
"VCUI", "VCUI", "VCRA", "VCRA", "VCVD", "VCVD")

Any suggestion would be helpful

Without regex, split on "_", drop last item, then paste it back again:

sapply(strsplit(x, "_"), function(i) paste(head(i, -1), collapse = "_"))

Or with regex (stealing the code from tools::file_path_sans_ext function which is used to drop file extensions):

sub("([^_]+)\\_[[:alnum:]]+$", "\\1", x)

# [1] "BAF60A_KD_1"     "BAF60A_KD_2"     "BAF60A_RA_1"     "BAF60A_RA_2"     "BAF60A_VD3_1"    "BAF60A_VD3_2"
# [7] "BAF60C_KD_1_UI"  "BAF60C_KD_2_UI"  "BAF60C_KD_RA_1"  "BAF60C_KD_RA_2"  "BAF60C_KD_VD3_1" "BAF60C_KD_VD3_2"
# [13] "VC_1_UI"         "VC_2_UI"         "VC_RA1"          "VC_RA_2"         "VC_VD3_1"        "VC_VD3_2"


Snakemake: Migrating from deprecated cluster.json to new profiles.yaml I am an avid user of Snakemake. Recently we have been refreshing our pipelines and I saw that a cluster.json file is no longer the recommended way to store the cluster configuration.

I used to start my pipelines like this:

snakemake --cluster-config cluster.json --cluster "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}" --jobs 245 --use-conda --rerun-incomplete --keep-going --latency-wait 60

With the cluster.json looking like this:

{
    "__default__" :
    {
        "cluster"   : "qsub",
        "jobs"      : 70 ,
        "nodes"     : 1,
        "ppn"       : 1
    },
    "rsem_quantify" :
    {
        "ppn"       : 8
    },
    "star_map_human" :
    {
        "ppn"       : 4
    },
}

Now, I understand that some of this goes in to the new profiles.yaml like this:

cluster: qsub
jobs: 70

But how do I set the ppn values per rule? And how do I set the other configuration options into profile.yaml (like: snakemake "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}", etc.)

Any examples are appreciated, https://github.com/snakemake-profiles is not really helping

It looks like you can just bundle your cluster config in its own YAML file, and reference that from within a profile YAML.  Here's a clumsy little example using a fake job submission script: https://github.com/ressy/example-snakemake-profiles

The current documentation makes it sound like cluster configuration support is going away and being replaced by profiles:

  While still being possible, cluster configuration has been deprecated by the introduction of Profiles.

...but from those example profiles, that might be misleading; here's an example template for a profile YAML that references a separate cluster-config YAML:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/config.yaml

...
cluster-config: "$((INSTALDIR))/overwrite_cluster_param.yaml" #abs path
...

That cluster-config YAML has the same default and per-rule entries just your JSON version does:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/overwrite_cluster_param.yaml

...
# specific parameters for certain rules, which need more time/memory

#run_assembler:
#  queue: bigmem
#  time: 100
#   threads and memory definen in config file

So I think you can just set up a pair of YAML files accordingly including options for your rsem_quantify and star_map_human rules.  From my own test it looks like it handles the cluster argument templating just as it did before when using a profile.  (If this is correct we should nag them to make the new documentation more explicit though.)


What is the index fastq file (sample_I*.fastq.gz) generated when demultiplexing Illumina paired-end runs? What is the index fastq file that comes with some Illumina sequencing datasets? (The samplename_I*.fastq.gz file.)

For example, I recently received some 10X Chromium reads for two libraries sequenced on the same lane. This was a 2x150 sequencing run, so there should be two fastq files. The sequencing center demultiplexed the libraries and generated two separate directories - one for each library.

In each directory there are three fastq files:

Mysample_I1_001.fastq.gz
Mysample_R1_001.fastq.gz
Mysample_R2_001.fastq.gz

I know that the R1 fastq file contains the forward reads and the R2 fastq file contains the reverse reads, but what is the index file? The first few sequences do not match the 10X indices for this library (ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA).

>bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | head
NTTGGGTA
NGCAATGG
NAGCCCAC
NAGCCCAC
NCATTACT
NCATTACT
NGCAATGG
NAGCCCAC
NTTGGGTA
NAGCCCAC

tldr - The I*.fastq.gz file contains the read index sequences.
long explanation
Illumina uses a program called bcl2fastq to demultiplex sequencing runs.
This software takes a list of samples and their associated indices and uses those sequences to make one or more fastq files per sample, binned by one or two index sequences on either end of the sequencing molecule (i5 and i7 indices, see page 6 for HiSeq).
Illumina sequencing isn't perfect though, and there are sometimes errors in reading the index sequence. For example this index, CAGCCCAC can easily have read errors in the A sandwiched between many Cs: CAGCCCAC -> CAGCCCCC.
Instead of throwing out all of the reads with indexing sequencing errors, the bcl2fastq program includes reads that are discernibly derived from a true index for that sample as long as there is no overlap with another sample.
In the example provided above, the four indices for mysample were: ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA. If we look at all of the index sequences in the Mysample_I1_001.fastq.gz file, we will see that sequences with the four correct indices are indeed the most prevalent, however there are also reads that have indices with sequencing errors derived from the correct indices.
bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | \
    sort | uniq -c | sort -k1 -nr | head

count      index       source
41362311   CAGCCCAC    True index
37209190   GGCAATGG    True index
36863213   ACATTACT    True index
33674467   TTTGGGTA    True index
 1140358   NAGCCCAC    CAGCCCAC
 1026099   NGCAATGG    GGCAATGG
 1016754   NCATTACT    ACATTACT
  933342   NTTGGGTA    TTTGGGTA
  119626   TTTGGGGA    TTTGGGTA
   98657   GTTGGGTA    TTTGGGTA
   96625   GGCAATGA    GGCAATGG


Difference between isoforms and paralogs in transcriptoms? I've assembled RNAseq data into a transcriptoms using Trinity. There's a option to keep only the longest isoforms for each transcripts and it lead me to wonder how it deals with duplicated genes (-->transcripts). I haven't found an anwser in the original paper. Althought I reckon paralogs and isoforms are two distinct features I don't see a way to make a distinction when you only have trasncripts.

If a reference genome assembly exists, then reference-guided assembly with Trinity will keep paralogs, because it treats each genomic region as a different assembly set.
If a reference genome doesn't exist, paralogs will probably result in the construction of non-representative "isoforms" containing variants present in multiple genes. Inchworm does a greedy assembly of reads, so if a chimeric transcript can be constructed from the reads (which is likely with paralogous genes), then it might end up being fed to the subsequent steps. Whether that is actually emitted at the end by Butterfly will mostly depend on how the paths through the putative assemblies are traversed, which probably depends on the tolerance of Bowtie to errors.
If paralogous transcripts are sufficiently similar that the differences are within the error tolerance for assembly, then I would expect that they'd be collapsed into a single transcript.


PCA plot did not work in single cell RNA-seq I run plotPCA for single cell RNA-seq data, while I get error message (I use R 4.0). I attached the code and error message here. Did anyone know the reason and how to deal with it. Thanks!

library(pcaMethods)
library(SC3)
library(scater)
library(SingleCellExperiment)
library(pheatmap)
library(mclust)
set.seed(1234567)
deng <- readRDS("deng-reads.rds")
plotPCA(deng, colour_by = "cell_type2")

Error in value[[3L]](cond) : invalid subscript 'type' in 'reducedDimNames(<SingleCellExperiment>)'
'PCA' not in 'reducedDimNames(<SingleCellExperiment>)'

I also get error message when I run the code as below:

colData(deng)$tSNE_kmeans <- as.character(kmeans(deng@reducedDims$TSNE, centers = 8)$clust)

# Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

The error is pretty obvious. There is no PCA data in the respective slot of the SingleCellExperiment that plotPCA could use, at least it is not called PCA which by default is expected. Use runPCA first, please read the manuals of scater and scran first, they have step-by-step code.


which NCBI tool is optimized to identify a species from a DNA fragment? “The Iceman” was a man who lived 5300 years ago and whose body was recovered from an Alpine glacier in 1991.  Some fungal material was recovered from his clothing and sequenced.

Ice man : found as a frozen corpse in the Val Senales Glacier in Italy.

Thought to be carrying Fomitopsis betulina (formerly known as Piptoporus betulinus) commonly known as birch polypore, birch bracket, or razor strop) is a common Basidiomycota brown rot macrofungus growing on decaying birch wood.

Other species known : Fomes fomentarius (commonly known as the tinder fungus, false tinder fungus, hoof fungus,] tinder conk, tinder polypore or ice man fungus)

Given two DNA fragments from NCBI genbank, try to find what species these DNA fragments belongs to or what species these fragments are most closely related to.

>Z54169.1 Iceman fungal clone 'Sim C-a1' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAAACTTGGGCCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTTCCTTCTGGAGAACCTCATGCCCTTCACTGGGTGTGTTGGGGAA
CCAGGACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAAT
AATAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCG
GGGGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCAT
TTGCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAG
TCTTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGA
GAAATCAAAGTTTTTGG

>Z54155.1 Iceman fungal clone 'Sim C-a2' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAACCTTGGGTCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTCCTTCTGGAGAACCTCATGTCTTTTATTAGGCGTGTTGGGGACC
AGAACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAATAA
TAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCGGG
GGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCATTT
GCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAGTC
TTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGAGA
AATCAAAGTTTTTGG

Now I know I can just google what species of mushrooms (which is where the above background info is from). But say that I was unable to and had to just go off the DNA fragments in order to find the species.

I've tried BLASTn, but I'm not sure if it is the right tool or if I need to properly calibrate the parameters to use it.

so far, I've tried:

Under “Choose Search Set:
Database : Non-redundant protein sequences (nr)
Organisms excluded: Animalia (taxid:33208), Plantae, Prokaryota  (taxid:2), Protozoa
Organisms included: fungi  (taxid:4751)
Included: Models (XM/XP),   Uncultured/environmental sample sequences
Under Program Selection:  Optimized for Highly similar sequences (megablast)

Algorithm Parameters (General Parameters, Scoring Parameters, Filters and Masking,
PS/PHI?DELTA BLAST)

are all kept default.
```

My conclusion from the tree is the fungi appears to be genuine and not an extant contaminant (which is very easily done). It is not Fomitopsis betulina because this is division Basidiomycota. The fungi here is division Ascomycota and these are very distinct from Basidiomycota.

There have been many examples of contamination in ancient DNA work, but this doesn't appear to be one of them. The rationale is given below.

The first fragment was not interesting and described 100% identity to ascomycetes. The tree above is a simple distance tree from the second fragment of 18S which shows it is an outgroup to the ascomycetes (green terminal taxa). The ancient fungi contains 3 (I think) unique nucleotide residues that have not been sequenced in any extant ascomycetes to date and sequences are unique to any known species. The location of the 3 unique nucleotides is tightly clustered within the locus with 100% identities either side. This is significant because it demonstrates it is not a PCR mosaic artefact (yep its been done in ancient DNA work and subsequently proven artifactual). In the latter PCR mosaic artefact I'm referring to (which I wouldn't cite), was a fungi/bacterial mosaic from memory.

Given the homology of 18S in the extant ascomycetes I would say this appears to be a new genus, however you would need to be sure this was not sequencing error and this is very old DNA, which is subject to degradation even when frozen. The implication of the tree is quite unusual, it implies widespread extinction of this group of fungi in recent history (5000 years is slight in evolutionary terms) BUT it could be a sampling bias, i.e. no one has sequenced those fungi and NCBI is not representative.

I assume that you have performed phenotypic/morphological analysis on the spore to believe it is Basidiomycota.

The fragment in yellow is the query and the sequence next to it was the sequence already deposited on Genbank of this ancient isolate. You would need to look up the biology and ecology of ascomycete to understand its significance, but it is in a lineage of its own, so maybe difficult to interpret by comparison to extant species. You could concatenate the two alignments based on each fragment, but the problem with 18S is that is a hard to model and you'd need to perform secondary RNA structure modelling to assess whether these mutations would be viable.

In conclusion, it definitely looks genuine but it is not the species you think it is, in fact not even close.


ENA minion 2D data doesn't actually have 2D data I am looking at some yeast ONT data from De novo yeast genome assemblies from MinION, PacBio and MiSeq platforms, Giordano et al 2017.

I noticed that the Minion fastq files have both a forward and reverse file. The methods paper says that these reads are from the SQK-MAP006 and SQK-MAP007 protocols which both appear to be 2D. Are these just the 2D forward and reverse files?

Still being unsure, I checked if this was true by counting the number of characters in each read.

>> zcat ERR1883385_1.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/1

>> zcat ERR1883385_2.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/2

>> zcat ERR1883385_1.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1675

>> zcat ERR1883385_2.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1693

As you can see, the read is the same but the forward and reverse strands are different lengths.

Is there some way to just download the poretools equivalent to  2D from ENA without having to download the fast5 files and extract them myself?

Probably not. ENA doesn't have any ability to specifically flag a file as 2D, and the nanopore file formats have changed a lot, so you pretty much need to rely on the submitter to separate out the reads in a meaningful way.

Also, there wasn't really any consistent convention for what ONT used in their read headers. Looking at some of my old files they had a plain sequence ID (i.e. without anything suggesting temp/comp) in different template/complement folders in the FAST5 files for the components of a 2D read and the consensus sequence, but later changed it to _2D_template and _2D_complement as in the example you have (possibly with an intermediate _1D_template and _1D_complement format). Parsing/understanding those FAST5 files can be difficult, which is why poretools was created in the first place.

It's still a bit of a mess, with the 1D$^2$ reads ending up in a Frankenstein-like state that's partially 1D and partially 2D.


How should the popular press compare similarity of genomes? Note this is a question from a lay reader.

I've read in the popular press that ~1-4% of the genome of non-African Homo sapiens is inherited from Neanderthals, or that Melanesians derive a similar amount from the Denisovans.

This seems like a lot, until you read (elsewhere!) that Neanderthals are 99.7% identical to humans, so most of the fraction inherited is identical in any case.

Or taking it another way, reading that two individuals will typically have 99.9% identical DNA seems pretty similar, until you realise you are only 0.3% 'different' from a Neanderthal, or 15% 'different' from your dog.

It seems to me this use of percentages is confusing, misleading, and open to abuse. When comparing species (for example) I would think that they should speak of differences in terms of  base pairs that we could not otherwise have had.

How should the lay (or popular science) press best communicate similarities or differences in genomes when talking about individuals, races, or species?

1–4% is from an evolution point of view. 99.7% is from a sequence point of view. Because they are measuring different things, they can be compatible with each other. The correct interpretation is: 1–4% of a non-African genome is inherited from Neanderthals and the sequence of this 1–4% differs from modern human sequence by 0.3%. PS: so, for each non-African, Neanderthals contributed 1–4% * 0.3% of base pairs not present in modern human. This sounds tiny, but even a single SNP may lead to severe disease. This tiny fraction has known functional impact to diabetes and high-altitude adaption.

Others pointed out two related questions:

How many genes do we share with our mother?
Do apes and humans share 99% of DNA or 99% of genes? What is the difference?

These questions are really about the ambiguity in the definition of "gene". Importantly, traditionally in popular science and population genetics, a gene could mean an allele. This makes differences.

For example, by common uses of "gene" in biology, the answer to "how many genes do we share with our mother" would be nearly 100% – every human has nearly the same set of genes. However, clearly this is not what the asker wants. By gene, he/she really means allele. The expected answer is 50% – we inherit 50% of alleles from our mother.

In the end, what is a "gene" exactly? I don't know. Biologists can often figure out the meaning of word "gene" in the right context, but non-biologists will sometimes have troubles with the ambiguity.


Software recommendation: find DNA sequence distribution over entire transcript I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms.

I have already found the number of occurrences of the sequences I am interested in using the seqinr and biostrings packages. I have been using biomaRt to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together.

The issue I am having is that biomaRt returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?

Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.

Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.

library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>%
    unlist() %>%
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>%
    tolower() %>%
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"),
            filters = "mgi_symbol",
            values="Grin1",
            mart = ensembl,
            bmHeader = TRUE)

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns


Intersecting two different files with one "master" file based on different columns I have the following sets of data:

file1:

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

file2:

1   15776231    2   5008    G:5002  A:6
1   15776239    3   5008    C:3358  A:14    G:1636

file3:

1   16855677    2   5008    A:5003  C:5
1   16855689    3   5008    G:4957  A:41    T:10

Most lines have 6 columns but a few lines also have 7 columns which have data that I want to preserve.

Can a script be written which does the following:

if $2 of file2 comes in between $2 and $3 of file1 (between the interval of $2 and $3) and $2 of file3 comes in between $5 and $6 of file1 (between the interval of $5 and $6), then the output would be $1,$2,$3,$4 from file1 and $2,$3,$4,$5,$6,$7 from file2 and $5,$6,$7 from file1 and $2,$3,$4,$5,$6,$7 from file 3, all in one line. So something like this:

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776239    3   5008    C:3358  A:14    G:1636  16855676    16855696    CTGTCCAGCAGAGGGCGGTG    1   16855689    3   5008    G:4957  A:41    T:10

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2 and file3. The intersection is based on $2  in both the files, and all the other columns after it are just copied

(There could be many lines from file2 and file3 intersecting in the same line in file1,so that line could be repeated based on the number of intersections)


tximport txidcol and geneidcol selecting wrong column I'm trying to run a differential expression analysis of 2 samples of the same species, but with different treatments using DESeq2. I'm using Ensembldb to create a tx2gene dataframe for tximport, however tximport does not select the columns I specify, and instead always selects the tx_biotype leading to the following output:

.                    25_quant    28_quant
nontranslating_CDS    10422.41    10022.74
protein_coding     22196722.59 24002269.26

I'm running tximport using the following params:

txi_tx <- tximport(files, type = "salmon", tx2gene = tx2gene,
                   txIdCol = "tx_id", geneIdCol = "gene_id")

I've gotten around this by deleting all columns from my tx2gene DF apart from tx_id and gene_id. However, I think this might be causing problems in the DESeq2 analysis somehow, and I would like to rule out that possibility.

This is what my names(tx2gene) looks like:

"tx_id"            "tx_biotype"       "tx_seq_start"     "tx_seq_end"
"tx_cds_seq_start" "tx_cds_seq_end"
"gene_id"          "tx_name"

What am I doing wrong?

By setting type="salmon" you're overriding your settings for things like txIdCol. The various type options (except none) are presets to set things correctly for the output from various programs like salmon or kallisto. If you need to specify custom column names, set type="none".

As an aside, if I'm interpreting your question correctly you have only two samples in total. You can not perform a meaningful analysis with an unreplicated experiment, so if that is indeed your case then don't waste your time further with this dataset.


PDB format: remark number for free text I would like to add a text to PDB files that I'm processing with my tool, rna-pdb-tools. Someone points that the way I'm using it right now it's not correct (https://github.com/mmagnus/rna-pdb-tools/issues/48).

I use HEADER right now which is not correct.

HEADER Generated with rna-pdb-tools
HEADER ver 37c5b4e-dirty
HEADER https://github.com/mmagnus/rna-pdb-tools
HEADER Mon Oct 10 22:18:19 2016
ATOM      1  P     C A   1     -19.687  -3.296  65.469  1.00  0.00           P

Do you happen to know which remark number to use? (http://www.wwpdb.org/documentation/file-format-content/format33/remarks1.html)?

From this document, it looks like REMARK 250 is the way to go.


How to convert BED to GFF3 I would like to convert a BED format to GFF3.

The only useful tool that I could find via a google search seems to be Galaxy, and I do not feel very comfortable with online tools, plus the webserver is currenlty under maintenance.

Does anyone knows about a command-line tool that can handle this conversion?

Edit: here are some lines of my BED file:

$ head -4 last_minion-r7_sort.bed
211000022278137 175 211 8e5d0959-7cdb-49cf-9298-94ed3b2aedb5_Basecall_2D_000_2d 42  +
211000022279134 0   503 e8a9c6b8-bad2-4a7e-97d8-ca4acb34ff70_Basecall_2D_000_2d 69  -
211000022279134 24  353 e258783d-95a3-41f5-9ad5-bb12311dbaf4_Basecall_2D_000_2d 45  -
211000022279134 114 429 26601afb-581a-41df-b42b-b366148ea06f_Basecall_2D_000_2d 100 -

The bed file thus has 6 columns as for now: chromosome, start coordinate, end coordinate, read name, score, strand. This file was obtained from conversion of MAF format (as output of alignment of RNA-seq reads to reference genome, using LAST) converted to SAM using maf-convert, then to BAM using samtools, finally to BED using bedtools.

The aim of my conversion is basically to convert SAM -> GTF, for post-processing. Since there is no straightforward way to do this, I am going through steps, the only way to do this in my knowledge is : SAM -> BAM -> BED -> GFF3 -> GTF but for now I am stuck in the BED -> GFF3 part.

Galaxy has API and API-consuming libraries (such as BioBlend) that will allow you to interactively script against it without opening the graphical interface at all.

However you can also take almost any tool out of Galaxy and use it independently since everything is open source. The converter you mentioned is available as a Python script here and the tool 'wrapper' which you can use to understand how to invoke the python script is next to it.


about GO terms's name Did anyone know whether the GO terms can include more detail information? Like I can get the DotPlot of the GO terms as below:

The problem is that some of the genes, in the GO terms, is more about sinoatrial-node (SAN). So whether there is any method that we can do to dig deeper into the GO terms?

The various Gene Ontologies are arranged as Direct Acyclical Graphs. This means that each term can have children and parents. When a gene is annotated with a GO term, it also inherits all of that term's parents. Therefore, when you run a GO enrichment analysis, it is very common that you find the broader terms enriched since those are the ones shared by most genes.

For example, the Biological Process (BP) GO term GO:0003163 sinoatrial node development has the following tree:

￼ GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

This means that any gene annotated with GO:0003163 will also be annotated with every one of the terms above it in the tree, all the way to the top of the graph. All genes with at least one BP GO term will also be annotated with GO:0008150, the root of the tree.

As you can see, your "muscle tissue development (GO:0060537)" is also there and is a parent of "sinoatrial node development". However, if you just ran an enrichment analysis (I ave no idea what you did since you haven't actually told us anything), then what will have happened is that there are many genes annotated with various different child terms of GO:0060537, but since all child terms will inherit the parent annotations, what you see is the parent and not the specific term you are looking for.

So yes, there are far more specific terms that what you show, but whether or not they appear in your results will depend on what you are doing.


How can I add for several bed files the header : track type=narrowPeak name=“narrowPeak” preferably in python ,can handle with R I want to create custom tracks from these files
I can add the line : track type=narrowPeak name=“narrowPeak” manually by opening it with text editor:
track type=narrowPeak name=“narrowPeak”

but I wonder if there is no simple way to do it with code- to open all files that end with bed and add them this line.
Didnt find this option in bed tools, If it is possible to add the name extracting it from the file name it would be amazing.
Many thanks

I don't think you need to go to R or Python to do that. In bash
set -euxo pipefail
for f in *.bed
do
  bn="$(basename ${f})"
  on="${bn%%.bed}_with_header.bed"
  echo 'track type=narrowPeak name="narrowPeak"' > "${on}"
  cat "${f}" >> "${on}"
done

EDIT: I just read you want to set the name based on the filename.
EDIT2: And also to add an arbitrary directory
set -euxo pipefail

INDIR="/path/to/bedfiles"

for f in $(find "${INDIR}" -name "*.bed")
do
  bn=$(basename "${f}")
  on="${bn%%.bed}_with_header.bed"
  dn=$(dirname "${f}")
  echo "track type=narrowPeak name=\"${bn}\"" > "${dn}"/"${on}"
  cat "${f}" >> "${dn}"/"${on}"
done


scRNA-seq differential transcript usage Many of the modern gene-quantification tools (Salmon/Kallisto) output transcript-level (as opposed to gene-level) data. All of the scRNA-seq analysis I have seen just uses the gene-level values. I suppose I can just use the transcript-level matrix, but then I lose any links between the different isoforms of the same gene. Maybe it does not really matter?

There are a few options for bulk RNA-seq and there was a similar question here, but are there any single-cell options?

I was hoping to find some tool that already exists.

The short answer is, none exist yet. I see a number of hurdles to overcome first.1

For example, no one has yet characterized saturation for transcript discovery, i.e., the reads per cell needed to maximize transcript information recovery.  This is a typical stepping-stone as a sequencing approach matures, and until it's done, it makes it unlikely that most scRNA-seq datasets are going to be appropriate vis-à-vis statistical power, even if there was a tool.  That is, if people are only aiming for cell-type discrimination via gene-level markers, datasets will likely be underpowered to measure transcript level differences.

Another issue is the lack of any transcript-level scRNA-seq simulators. All the simulators out there that I know of produce gene-level matrices.  Having a transcript-level simulator is crucial for benchmarking and characterizing the performance of differential expression calling methodologies.  Without anything to test against, that leaves a huge blindspot in terms of estimating the uncertainty (false positive rates, sensitivity) that any tool might have.  However, since some scRNA-seq models (e.g., Velten, et al. 2015, or Lopez, et al., 2018) are generative, perhaps this is shouldn't be too difficult to adapt to generating FASTQs.

Simulators are fine and all, but ground truth is probably the key hurdle of the enterprise right now. While there are definitely groups using smFISH to distinguish isoforms quantitatively on a single-cell level, there doesn't appear to be enough smFISH data out there to measure a transcriptome-wide scRNA-seq method against. Hopefully this will change in the coming year or two, as results from more highly-multiplexed single-molecule technologies become available.

1 This is just my opinion formed as I try to keep up with the literature, so I could be missing things.  Comments/corrections certainly welcome!


How to calculate Gene Ontology terms in python I am testing different Protein-Protein Interaction networks alignment tools. I observed that all alignment tools show some Gene Ontology terms similarity. How can I calculate Functional similarity / average functional similarity for Gene ontology's terms Biological processes, Cellular components, Molecular functionality of any two proteins networks

You can calculate gene ontology similarity with the GOSemSim R package (paper)(disclaimer I'm a contributor to this package). You have several similarity scores implemented, some of them are across the three subontologies.

In python you have the GOATOOLS (paper). I think there was another Python package but now I don't remember its name.


REMC SNP annotations for transcriptome prediction - epiXcan I have been trying to use a pipeline for predicting gene expression for a target gene (Although you can do it for many) - the name of this method is epiXcan. Here is a link to the paper:

https://www.nature.com/articles/s41467-019-11874-7#Sec22

One of the very first steps in the pipeline is obtaining SNP 'annotations' from an epigenome reference consortium known as the REMC. In the methods paper they say:

"We first prepare eQTL statistics (computed with MatrixEQTL4) and SNP annotations (extracted from REMC https://egg2.wustl.edu/roadmap/web_portal/). For each eQTL tissue, we use the matched REMC tissue to
extract the corresponding annotations".

I have completed the first step of obtaining eQTL statistics. However obtaining SNP annotations is where I am lost.

In my case, I am only trying to predict the expression of one gene (DRD4) and so my region of interest is the pre-frontal cortex. If you visit the website they provide though, it doesn't seem as though annotating a large list of SNPs - at least in the way they are presenting - is possible.

Can anyone provide some insight or direction on how to proceed? How does one obtain SNP annotations from the REMC? I have tried emailing the authors of the paper and they have not been responsive.

Thank you,

After much research I found out how to obtain the annotations.

You must go to the REMC website: https://egg2.wustl.edu/roadmap/web_portal/

Go to the chromatin-state-learning tab, and choose the n-chromatin state model desired - the core has 15 states, then 18, then 25 (imputed states).

Following this, you must go to the download URL link and, you will find a lot of files, primarily with the following coding information: E0*.. This corresponds to a particular tissue/cell type and so you must choose the one relevant to your research question. The way they are correspond to each tissue/cell type can be found in the following paper:

https://www.nature.com/articles/s41467-019-11874-7#Sec22


Where can I find Single Cell Data with Location "Coordinates"? Does single cell data typically have the following meta-data: the "coordinates" (e.g. on a tissue, adjacent tissues) saying where each cell in the sample was located relative to other cells? If not, is it possible to reconstruct this with other meta-data on the cells?
With the ultimate goal of working hands on with such location-tagged data, I am hoping for suggestions on the correct terms to search for, and even references to previous studies and/or easy to use public datasets.
For visual reference of the idea I have in mind, consider this picture:

I expect expression of any given gene in the skin cells in Condition 1 to be very different according to the location of the cell. For condition 2, I wouldn't expect the expression of a gene in the skin cells to be overly different, if at all. I want to see if I can formalize this idea using data mining techniques that were created for other purposes, but first I need the proper data.
edit: Here is an example of the format I eventually hope to work with.
Fake table 1: Each row represents measurements taken for a unique cell, while the last three columns have coordinates for its location.
\begin{array}{r|lllllllll}
  \hline
 cell & phenotype1 & phenotype2 & \dots & gene1.expr & gene2.expr & \dots & x.loc.coord & y.loc.coord & z.loc.coord \\
  \hline
1 & 0.69 & 1.34 & \dots & 1.91 & 0.21 & \dots & 1.12 & 0.05 & 1.09 \\    2 & 0.34 & 0.92 & \dots & 1.74 & 2.03 & \dots & 0.57 & 0.46 & 0.24 \\
   \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
  n & 1.97 & 1.3 & 0.96 & 0.19 & 0.66 & \dots & 0.25 & 0.02 & 1.27 \\
   \hline
\end{array}
Alternatively, if the (x,y,z) coordinates of Fake Table 1 are unavailable or not reconstructable, are there datasets which help us construct an adjacency list of cell pairs according to their location, such as Fake Table 2?
Fake Table 2: An adjacency table with pairs of cells which were next to each other when the measurement was taken.
\begin{array}{r|cc}
  \hline
pair & cell1 & cell2 \\
  \hline
1 & 1 & 2 \\
  2 & 1 & 3 \\
  3 & 1 & 4 \\
  4 & 2 & 3 \\
  5 & 2 & 4 \\
  \vdots & \vdots & \vdots \\
  m-1 & n & n-2 \\
  m & n & n-1 \\
   \hline
\end{array}
I come from a statistics background with a basic understanding of this type of data from talks I've attended, but not hands on experience. I mostly just want to explore this type of data to inform my future work. I am open to any type of scRNA-seq data with cell location "coordinates," if that's the correct thing to ask.

Does single cell data typically have the following meta-data: the
"coordinates" (e.g. on a tissue, adjacent tissues) saying where each
cell in the sample was located relative to other cells?

No. Typical scRNA-seq is just capturing random cells in a tube with no additional information. The technology you are looking for is spatial transcriptomics where you are measuring RNA levels at a particular location, but that is still not on the single-cell level.

If not, is it possible to reconstruct this with other meta-data on the
cells?

There have been efforts to reconstruct spatial relationship. For example: Satija et al. 2015. There should be more recent approaches that I am not aware of.


Comparing phylogeny in R So I want to compare the phylogeny created using two methods for example Maximum likelihood and maximum parsimony.Is there any way to compare the two phylogeny ?

I did read about phangorn but not sure if its the right R library for comparative analysis.

Any suggestion or helped would be highly appreciated

My data file

library(phangorn)
library(phytools)
library(dendextend)

data <- read.dna("abhi_seq/clean_dup_align_fast.fas", format = "fasta")
data
dat <- as.phyDat(data)

dm <- dist.ml(dat)
treeUPGMA <- upgma(dm)
treeNJ <- NJ(dm)

layout(matrix(c(1,2), 2, 1), height=c(1,2))
par(mar = c(0,0,2,0)+ 0.1)
plot(treeUPGMA, main="UPGMA")
plot(treeNJ, "phylogram", main="NJ")

dev.off()

parsimony(treeUPGMA, dat)
parsimony(treeNJ,dat)

tr.mp <- optim.parsimony(treeNJ, dat)

#tr.ml = optim.pml(treeNJ, dat)

fit <- pml(treeNJ, dat)
fit <- optim.pml(fit, rearrangement="NNI")

fit.ini <- pml(treeNJ, dat)
fit.ini

fit <- optim.pml(fit.ini, optNni=TRUE, optBf=TRUE, optQ=TRUE, optGamma=TRUE)
fit

tr.ml <- root(fit$tree,1)

tr.mp.ultra<-force.ultrametric(tr.mp)
tr.ml.ultra<-force.ultrametric(tr.ml)

is.ultrametric(tr.mp.ultra)
is.binary.tree(tr.mp.ultra)
is.rooted(tr.mp.ultra)

dd.ml.ultra<-as.dendrogram(tr.mp.ultra)

  Error in ape::as.hclust.phylo(object) : the tree is not rooted

I ran into this error this error which says trees is not rooted

Comparing phylogeny code updated ** Working thanks to ***thomas duge de bernonville* thomas for putting codes together and fixing the errors

library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

a<-read.alignment("abhi_seq/clean_dup_align_fast.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a.rooted<-root(pars.a, outgroup="AAA64460", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))

mt <- modelTest(a.phydat, tree=upgma.a,multicore = TRUE,mc.cores=10)

#ml.a = pml(upgma.a,a.phydat)
#fitJC <- optim.pml(ml.a, TRUE)
#ml.a.rooted<-root(midpoint(fitJC$tree), outgroup="AAA64460", resolve.root=T)
#ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

#########################################################################

######################################################################33

#mt <- modelTest(dat, tree=tree, multicore=TRUE)
mt[order(mt$AICc),]# choose best model from the table according to AICc
bestmodel <- mt$Model[which.min(mt$AICc)]
env = attr(mt, "env")
fitStart = eval(get(bestmodel, env), env)
fit = optim.pml(fitStart, rearrangement = "stochastic",optGamma=TRUE, optInv=TRUE, model="GTR")#tree
bs=bootstrap.pml(fit, bs=25, optNni=TRUE, multicore=TRUE)

ml.a.rooted<-root(midpoint(fit$tree), outgroup="AAA64460", resolve.root=T)
ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

Another 2 tree compare tools: ggtree and blatic 3


Proteins with one SS bond? I would like to find proteins with exactly one SS bond. Is there a database where I can search this? I've tried advanced search on https://www.rcsb.org/, but no such option, at least I could not find it.

In cases when you cannot find answer through the web interface, you can download all the data and search it locally.

Relying on the PDB annotations of SS bonds, you could search for PDB files that have exactly one SSBOND record:

$ cd $PDB_DIR/structures/divided/pdb/
$ zgrep -c ^SSBOND aa/pdb*.ent.gz | grep :1$
aa/pdb2aaq.ent.gz:1
aa/pdb3aad.ent.gz:1
aa/pdb5aap.ent.gz:1
aa/pdb6aa4.ent.gz:1
aa/pdb6aa5.ent.gz:1

All PDB files (gzipped) take about 30GB of disk space.

The caveat is that some entries are available only in the mmCIF format. Using mmCIF would be slighty more complex: you would need to count _struct_conn.conn_type_id being disulf. And you would need to download 10GB more.


How to modify DNA evolution model to fit actual data? I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).
I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms.
I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.
The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?

Mutation is not a random process because there are many mutational biases, e.g. the transition vs. transversion ratio. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

INDELible http://abacus.gene.ucl.ac.uk/software/indelible/
Evolver https://www.drive5.com/evolver/ (not to be confused with the PAML package of the same name)


Batch detection of CRISP proteins in fasta file Probably a naive question. I am inexperienced.

I am interested in identifying potential CRISP (Cysteine-rich secretory proteins) in a certain tissue transcriptome (ca. 20k sequences in fasta). I have detected signalP and estimated % of cysteine in sequences. However CRISP rely on a certain pattern, and I am not aware of batch-detected algorithms.

Please would anyone know the best/simplest way of detecting CRISP candidates?

[EDIT]

I have been asked for the pattern. I summarise what the literature says:

Several of them have a CAP domain;
There should be a region called 'cysteine-rich domain' (CRD) upon the carboxyl terminal half of the protein containing 10 of roughly 16 conserved cysteines.

I am particularly interested in venom-derived CRISPs, but I am not sure whether they cluster in a certain CRISP family.

References:
Roberts et al. (2007) Structure and function of epididymal protein cysteine‐rich secretory protein‐1.
Yamazaki & Morita (2004) Structure and function of snake venom cysteine-rich secretory proteins.

Is your goal to detect the actual proteins from the sequences? If you'll settle for detection of transcripts from the genes encoding these proteins as opposed to detection of the proteins themselves, here is an approach.

Aligning your transcriptome to a reference genome (probably with STAR) will be the first step. You have lots of short reads in your FASTA, but in their current state they are no more informative than a box of puzzle pieces. Aligning your reads to the reference will allow you to extrapolate which genes that you short reads correspond with. From here you can quantify the features (genes) in your transcriptome.

Lucky for you, there are only 3 genes in the cysteine rich secretory protein family, CRISP1, CRISP2, and CRISP3. Alignment with a reference (and annotations) and quantification of genes that your reads align to will give you all you need. You'll just have to find the three genes in your final counts table.

I recommend STAR as an aligner with the --quantMode GeneCounts argument. It's very quick if you've got the resources and is geared towards RNA-seq reads.


RepeatMasker annotation I collected my own denovo library using many different tools Pipeline (structure-based ) and I updated all the headers corresponding to RepeatMasker format, to obtain Repeat coverage in the genome and .tbl annotation.

when I run repeatmasker the annotation results assume to be Homosapien, where Homo annotation table ( .tbl) is limited as shown below Table (1)

Actually my denovo repeat library consists of more than what is in the table above, example, rolling circles, Copia, Gypsy etc..

RepeatMasker -Lib myOwnLibrary.fa Mygenome.fa

Table (1)

My question is how can I run RepeatMasker with my own library using -lib Flag and result like the table (2) below,

the table below assume to be vertebrates genome, this is my aim.

RepeatMasker -Species Vertebrates Mygenome.fa

Table (2)

Regards

RepeatMasker comes with a buildSummary.pl script.

You can input the .out file from your RepeatMasker run and there is a flag to enter your species

e.g.

buildSummary.pl -species Vertebrates rmrun.out > mysummary.tbl

In the mysummary.tbl, the summary will not assume your run to be Homo (default).


A good tool for gene locus visualization Does anyone know of a good tool for visualizing gene loci? Basically, I'd like to show how a group of genes look in different strains of a bacteria. I have annotation of course, I just need a tool to plug it into and get an image of the locus to compare.

I'm talking something basic like this:

I think in R, Gviz is the best choice. I have used it for human sequences but it will probably also be suitable for Prokaryotes.

Another option is MS powerpoint, if you only need the arrows...


snakemake - using output file not explicitly declared as input for subsequent rule I am writing a snakemake that, up to now should

use raw sequencing reads to produce a draft assembly using canu
use minimap2 to align the raw reads against the draft assembly, generated by canu.

The following two rules are described below:

rule canu:
    input:
        fqdir + "{sample}.fastq"
    output:
        directory("{sample}/canu")
    shell:
        "/home/bhinckel/downloads/canu-1.9/Linux-amd64/bin/canu genomeSize=4.1m -p {wildcards.sample} -d {output} -nanopore-raw {input}"

rule minimap:
    input:
        reads = fqdir + "{sample}.fastq",
        draft_assembly = "{sample}/canu/{sample}.contigs.fasta"
    output:
        "{sample}/aligned/{sample}.sam"
    conda:
        "envs/general.yaml"
    shell:
        "minimap2 -ax map-ont {input.draft_assembly} {input.reads} > {output}"

However, if I do snakemake -np I get Missing input files for rule minimap. Which does make sense, after all I am not specifying the input for rule minimap ({sample}/canu/{sample}.contigs.fasta) anywhere i.e. snakemake cannot know that the file {sample}/canu/{sample}.contigs.fasta will be generated by canu.

The thing is that canu produces several files as output (amongst which the file {sample}/canu/{sample}.contigs.fasta is one of them), that's why I specified its output as a directory.

One way of solving this I guess would be to specify all output files produced by canu as an individual rule output (e.g. {sample}/canu/{sample}.contigs.fasta, {sample}/canu/{sample}.report, {sample}/canu/{sample}.correctedReads.fasta.gz), though I am looking for a more robust way to do so (if that's possible at all).

Another way is to add one command after minimap: touch done.txt. This creates an empty file after minimap is finished. This can then be used in output rule for snakemake.

But specifying all the output files may not be a bad idea if these outputs are used in subsequent rules.


Why does this human bam file only have one copy of each chromosome? As we know that in human DNA sequence, one copy of chromosome comes from mother's DNA and another copy comes from father's DNA so as to form two copies of each chromosome in human DNA. So, if we extract exome sequence from its DNA then each exome must also have two copies of each chromosome (one from mother and another from father). But in either fastq file or bam file in WES data, I always find only one copy of each chromosome.

Q1: Where is the other gene copy in the sequence or have I have missed something?

Q2: How can I check the ploidy and find if there are two copies of each chromosome in WES data? How can I do this if my WES data are mapped to the reference and stored in a bam file?

The maternal and paternal copies of a chromosome are called haplotypes. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states.

Response to Q1

Your question, in other words, is: Why do .bam files not differentiate between haplotypes?

Your question gets to a more fundamental kernel of how most people do "genomics" today. When most people assemble a genome or create a reference genome, they actually just produce a fasta file that contains a single sequence for each chromosome. This, of course, is not biologically accurate as there are two very distinct sequences per chromosome. Most genomes only report a chimeric consensus between the two and call this a reference. These reference genomes are haploid, or haplotype-collapsed.

This is where your question comes in: You have mapped transcripts that biologically are derived from two different haplotypes to a reference genome that only contains one consensus sequence that (inadequately) represents both haplotypes. As a result, at a single locus in the bam file there will be mapped reads from both haplotypes. If your reference fasta file contained the exact sequence of both haplotypes, if you mapped reads to the reference, and if you only look at primary alignments, the reads will mostly map to their haplotype-of-origin.

This gets into another topic called phasing, wherein the order and orientation of most of the polymorphisms unique to each haplotype are determined using sequencing data. There are some problems with this as it relies on properly detecting variant sites. Software like GATK and others can find single nucleotide polymorphisms (SNPs) if there is good sequencing coverage, however detecting insertions and deletions is much more difficult. This gives a very SNP-biased view of the haplotype differences in any given genome. After finding variant sites, phasing software like  hapcut2 determines which variants fall on which haplotype and outputs blocks of variants that are thought to belong to the same haplotype.

Phasing alone is not enough to accurately reconstruct the exact sequence of both haplotypes due to the inability to detect all variants with read mapping. The gold standard of the future is diploid de novo genome assembly, in which each haplotype is assembled independently. This is an active area of research for people who develop genome assemblers and is tightly linked to the advances in PacBio and Oxford Nanopore long reads. This paper on trio canu is a good start for learning about one successful dioploid assembly technique.

Response to Q2

If you want to check ploidy of all of the chromosomes then you need at least 10X coverage of whole-genome shotgun data and a tool like smudgeplot and genomescope.

If you are trying to check local duplications or whole-chromosome duplications for something like a cancer sample you would also need whole-genome shotgun data. WES data does not give reliable information on ploidy since the read count at a given locus is dependent on the transcription level, not the actual quantity of homologous chromosomal DNA for the region in question.

If you are trying to (use a bam file to) look for evidence of chromosomal duplication or duplication of a locus in cancer samples you would look for an increase in whole genome shotgun coverage at that locus compared to a known normal locus. For example if all chromosomes but one have a coverage of 22, and one chromosome has an average coverage of 33, this is evidence of a trisomy. This logic can be applied to smaller regions as well (barring repetitive regions, paralogs, et cetera).


changing color key range to specified range in heatmap.2 function I have a tab separated text file as shown below

a   0.311   0.510   0.123   0.002   0.001   0.417   0.572   0.074   0.169
b   0.324   0.592   0.070   0.028   0.028   0.535   0.535   0.127   0.113

I am trying to use heatmap.2 function from gplots package and want to change the color key so that it ranges from 0 to 1

Following is my code:

d <- read.table("file.txt",sep="\t",header=F,row.names = 1)
colnames(d) <- c("-3","-2","-1","+1","+2","+3","+4","+5","+6")
d_matrix <- data.matrix(d)
colors<-seq(0,1,length.out=100)
heatmap.2(d_matrix,scale="none",col = redblue(75),trace="none",Colv = FALSE,Rowv = FALSE,breaks=colors)

But this gives me error:

Error in image.default(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 +  :
  must have one more break than colour

I want to change to color key such that it ranges from 0 to 1. How can I do it in the function heatmap.2

I added breaks argument to the heatmap.2 function as suggested from this post here but no luck..

You just need to add more colors. The following, for example, will work:

heatmap.2(d_matrix, scale="none", trace="none", Rowv=F, Colv=F, breaks=seq(0,1,0.01), dendrogram="none", col=redblue(100))

Note that length(seq(0, 1, 0.01)) is 101, therefore redblue() must return one fewer (100) colors (or colours if you prefer extra vowels).


How to calculate Centimorgan lengths from run lengths in one-to-one dna profile comparison For easier management of a database of dna profiles, I'm creating a tool to perform one-to-one matching on two profiles with snp genotype information.

Basically, something that works similar to GEDMatch's one-to-one matching feature.

Determining the length of runs of half or fully matched snp's, including some allowance for a genotyping error rate, is the easy part.

My question is about converting an snp run length to centimorgans.

I figure that the length needs to be adjusted by the likelihood of a change in the particular region of the chromosome.

Since my dna profiles and the positions of snps are based on Genome Reference Consortium Human Build 37, I figure that his information could be in the reference build, directly or derivable in some way.

Looking at the GRCh37 download files and the definitions of their formats, I'm unable to find such information.

Is this information present in GRCh37, or can it be downloaded from another reference source, and from which field in the dataset can it be retrieved or derived?

Try having a look at this link:

http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh37.map.zip

If you unzip the file, there is a file for each chromosome (1:22) and then 3 for the x chromosome. The first column of each file details which chromosome the file corresponds to, the 3rd column gives you the cumulative genetic distance in cM and the 4th column gives you the genomic coordinate (position) in bp.

You can then calculate the distance in cM between two positions from this data by subtracting the cM positions between the two positions you are interested in.


SRA Toolkit execution problem I am trying to install SRA Toolkit on MacOS, following all the instructions from NCBI website. It is installed but fastq-dump cannot be executed because "cannot execute binary files". I tried "chmod +x", noah from home-brew to make binary executable on MacOS but fast-dump still cannot be executed. I attached the screenshot of my terminal window. Can anyone help?

You can't run Linux binaries on OSX (or vice versa). Download Mac OSX binaries.


Does NCBI's blast API block my IP? I am trying to run a blatn command from both my laptop and within a google colab notebook. I am not sure why, but this command runs properly once and then, on the second execution, it gets stuck for hours without a reply.
Is it possible that NCBI blocks my IP? I've tried also to run web_blast.pl and this didn't work either.
biopython:
from Bio.Blast import NCBIWWW

fasta = '>seq0\n'+'AACCACCTGACGGGCTCTGT'
blast_res = NCBIWWW.qblast('blastn', 'nt', fasta)
print(blast_res.read())

web_blast.pl command:
web_blast.pl blastn nt input.fa > blast.out
Note, if this is due to throttling by NCBI, I am not sure how do they expect people to query their web service if they block you after one call.

You should check the User Guide:

The NCBI BLAST servers are a shared resource. We give priority to
interactive users. In order to ensure availability of the service to
the entire community, we may limit searches for some high volume
users. Interactive users of the NCBI webpages through a web browser
should not encounter problems. We will move searches of users who
submit more than 100 searches in a 24 hour period to a slower queue,
or, in extreme cases, will block the requests. To avoid problems, API
users should comply with the following guidelines:

Do not contact the server more often than once every 10 seconds.

Do not poll for any single RID more often than once a minute.

Use the URL parameter email and tool, so that the NCBI can contact you
if there is a problem.

Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays
if more than 50 searches will be submitted.


How to create this function for lapply? I have two lists of dataframes. List x and list my_list. Each having four dataframes and their names are same for the two lists, i.e., bray-curtis, Chao1, FB, and Shannon.
The dataframes within the two lists have same element names, but they contain different information. This is because list x has been generated after some operations on list my_list.
Now I want to write a function and run that on list x with `lapply. the command I want to write as a function is
forest(x$Chao1, slab = paste(my_list$Chao1$StudyID))

Here, the first Chao1 is coming from list x while the second  Chao1 is coming from the list my_list. It works for a single dataframe but I am unable to write the codes with lapply so that it runs across all the dataframes of the lists.
The desired outputs are four plots from four dataframes of list x.
Can anyone please help me?
Update: @gringer's approach works, but it saves all the plots in a single list. When I tried to save the images with following command it just saved the last image:
tiff(filename = "~/Desktop/image.tiff", width = 30, height = 30, units = "cm", res = 1200)

Is there a way to save all them individually?
Thanks

Assuming the lists are named, and you want to fetch the element with the same name from each list, I'm guessing something like this would work (or at least get you closer to a solution):
v <- lapply(names(x), function(var){
   forest(x[[var]], slab=(my_list[[var]])$StudyID)
})

When creating images that don't have a multi-page format (i.e. anything except PDF), and not changing the image name with each save, you need to add in a format specifier to tell R what to do with new created images. For example, this will add a three digit zero-padded number for each new TIFF image:
tiff(filename = "~/Desktop/image_%03d.tiff", width = 30, height = 30, units = "cm", res = 1200)


Build tree to calculate beta diversity I am trying to calculate beta diversity using the UniFrac distance on my data, each sample having OTUs values.
To calculate it, I saw it required building a tree with percentages that represented the OTUs.
However, I was unable to figure what percentage fit each OTU to build the tree from the instructions I found on the internet.
For example, I saw in this link: http://readiab.org/book/0.1.1/3/1#4.1.2,  that they build this tree:

Could you please explain to me how to build the tree?

Can you give a little more information about what tools you are using? How did you generate the OTU table etc?
FastTree is a very simple tool to generate a decent phylogeny. If you have the sequences of each OTU, you can align them in MegaX with muscle, or simply using muscle from the command line if you prefer. You can then generate the tree with this command (presuming you are on windows):
FastTree.exe -gtr -nt < aligned_sequences.fasta > otus_tree.tre

I hope that helps!


Read-informed statistical phasing I'm wondering if there are any tools out there that perform statistical phasing of genetic data (e.g. Eagle) that take read-based phasing from GATK into account.  I've found this paper which sounds like what I am looking for, but that tool was developed specifically for Hi-C data, and I only have exomes.

Thanks

shapeit4 is currently (probably) the most accurate and quickest phasing algorithm which can account for sequencing reads in the model. They suggest you extract the phase information from the reads using WhatsHap first however. The guide for WhatsHap can be found here and has some info on using data processed using GATK.


How to determine a protein's cellular location based on its sequence? I am wondering about the appropriate workflow to determine a protein's cellular location based on its sequence.

Let's say I have a sequence like this from a fasta file

MTIKEHKVVYEAHNVKALKAPQHFYNSQPGKGYVTDMQHYQEMYQQSINEPEKFFDKMAK
EYLHWDAPYTKVQSGSLNNGDVAWFLNGKLNASYNCVDRHAFANPDKPALIYEADDESDN
KIITFGELLRKVSQIAGVLKSWGVKKGDTVAIYLPMIPEAVIAMLAVARIGAIHSVVFAG
FSAGSLKDRVVDANSKVVITCDEGKRGGKTINTKKIVDEGLNGVDLVSRILVFQRTGTEG
IPMKAGRDYWWHEEAAKQRTYLPPVSCDAEDPLFLLYTSGSTGSPKGVVHTTGGYLLGAA
LTTRYVFDIHPEDVLFTAGDVGWITGHTYALYGPLTLGTASIIFESTPAYPDYGRYWRII
QRHKATHFYVAPTALRLIKRVGEAEIAKYDTSSLRVLGSVGEPISPDLWEWYHEKVGNKN
CVICDTMWQTESGSHLIAPLAGAVPTKPGSATVPFFGINACIIDPVTGVELEGNDVEGVL
AVKSPWPSMARSVWNHHDRYMDTYLKPYPGHYFTGDGAGRDHDGYYWIRGRVDDVVNVSG
HRLSTSEIEASISNHENVSEAAVVGIPDELTGQTVVAYVSLKDGYLQNNATEGDAEHITP
DNLRRELILQVRGEIGPFASPKTIILVRDLPRTRSGKIMRRVLRKVASNEAEQLGDLTTL
ANPEVVPAIISAVENQFFSQKKK

I can feed it to NCBI's blastp which then gives me a list of genes. I can check e.g. the first entry which also provides a link to SwissProt. There I finally see that this protein is located in the nucleus and also the cytosol.

My questions are:

1) Is there a more straightforward way of achieving this information?

2) Let's say someone annotated this protein to be a mitochondrial one. How could I check that i.e. (dis)prove that?

You could use Blast2go to find the predicted relevant gene ontologies of your sequences. If you look for the cellular compartment sub ontology you should find the location of such protein. I have never used this tool, but here is the paper describing how it works.

I wouldn't restrict to the first predicted match of blastp see if some more genes with a good fit are also located at the same sites. But otherwise your workflow seems right to me.

But to prove this is in X or Y you, or your team or collaborators, should always check with histochemistry if it really allocates on the predicted place. All the softwares predictions must be validated with some biological experiments to prove them.


Reading multiple raw files in Seurat I have multiple single cell samples to analyze and I'm following the instructions in Satija Lab's website. I want to merge all the count files from all the samples at once, and associate the metadata to each sample.

For the Read10x command, is there a way to read multiple files at once, such as merging all the filter_featire_bc_matrix files?

I have all the samples saved in a directory on my desktop, is it possible to call the directory within the script for it to read all the raw counts?

This is what I have thus far:

# load the required libraries
library(dplyr)
library(Seurat)
library(ggplot2)

-----------------------------------------------

# set up the working directory

wd = ("~/Seurat") setwd(wd)

# read all metadata

metadata = read.table("metadata_Jan2020.csv", sep=';', header=T, stringsAsFactors=F)
rownames(metadata) = metadata$LibraryID

#metadata
#mf = metadata

mf = metadata %>% filter(Type %in% c("PDAC", "NORMAL"))
mf = mf %>% filter(!(ID %in% c("PDAC-HTB2854")))`

Seurat has a vignette specifically for combining multiple 10x libraries.


bioawk removed part of FASTQ header I used

bioawk -cfastx 'length($seq) > 1 {print "@"$name"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

in order to keep a particular read length, but this command shortened the header  from

@A00199:161:HF3JLDMXX:1:1101:5882:1063 2:N:0:CTGTAGCC+GAACAGGC

to

@A00199:161:HF3JLDMXX:1:1101:5882:1063

How is it possible to fix the bioawk command to remain the full command?

bioawk -cfastx 'length($seq) > 1 {print "@"$name $comment"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

Note the addition of $comment, which is what follows the name on the first line of each entry.


Is there any publicly available multi-omics dataset? I am searching for a multi-omics dataset which may include genomics, transcriptomics and proteomics (eg. snyderome).
I need such a dataset for an introductory data-exploration purpose. So it can be from any host organism and any experimental methods/conditions.

Is there any database/repository where I can find it?

There is a database called OmicsDI, where one can search for multi-omics datasets.

Here's a link of the associated publication (Perez-Riverol, Yasset, et al. "Omics Discovery Index-Discovering and Linking Public Omics Datasets." bioRxiv (2016): 049205.) for more details.


Is it possible to obtain RefSeq transcripts from Ensembl database for GRCh37 to match Ensembl Genes? I know that with GRCh38, it is possible to extract the matching RefSeq transcripts for Ensembl genes/transcripts/proteins using Biomart from Ensembl.

However, for GRCh37, this option does not exists. Using the following MySQL command returns Ensembl transcript IDs for matching RefSeq transcripts:

mysql -u anonymous -h ensembldb.ensembl.org -A -P 3337 -D homo_sapiens_core_99_37 -e "SELECT transcript.stable_id, xref.display_label FROM transcript, object_xref, xref,external_db WHERE transcript.transcript_id = object_xref.ensembl_id AND object_xref.ensembl_object_type = 'Transcript' AND object_xref.xref_id = xref.xref_id AND xref.external_db_id = external_db.external_db_id AND external_db.db_name = 'RefSeq_mRNA';"

Also, I was trying to play around with the databases schema but I am not able to find any solution to problem.

Ultimately, I wish to obtain a table like the following:

Gene stable ID      Transcript stable ID    Protein stable ID   RefSeq match transcript
ENSG00000120329     ENST00000239451         ENSP00000239451     NM_031947.4
ENSG00000164219     ENST00000419445         ENSP00000404676     NM_005023.4
ENSG00000150594     ENST00000280155         ENSP00000280155     NM_000681.4
ENSG00000113578     ENST00000337706         ENSP00000338548     NM_000800.5

Thank you in advance.

The RefSeq match option in BioMart is from the Matched Annotation from NCBI and EBI (MANE) collaboration between RefSeq and Ensembl. It has only been calculated for the up-to-date gene annotation on GRCh38 so cannot be obtained on GRCh37. You can get mapping from Ensembl to RefSeq transcripts through BioMart as RefSeq mRNA ID (refseq_mrna in R) but this is not a perfect match like the MANE, it is a mapping based on sequence similarity and similar genomic location, and there can be mismatches between them.


Evaluating amateur bioinformatics results I've might have created a new algorithm for finding patterns in DNA. The technique has not been used before, and it managed to find a 33-mer with corresponding complement that exists both in E. Coli and in Salmonella. I think it's pretty neat but I don't know how to evaluate the results. In addition, I've developed this algorithm after a bioinformatics course where only the DNA sequences of E. Coli and Salmonella were provided, so I would like to test it on other sequences and have it evaluated by people who know what they are doing.

How would one get in contact with the right people for evaluating amateur bioinformatics result?

You need to understand the biological significance of the 33-mer locus.

I'd blast your 33-mer on NCBI and see what the biological significane is of the locus. If might be a transposon/ promoter region.
Look at the translation of the protein and if its part of a protein the location within the protein structure.
Blast the entire locus including 5' and 3' regions of the 33-mer target
Carefully assess the results for chance, e.g. conserved protein and the synonymous mutations just happened to drop out (Poison distribution).

Once you gain a greater understanding, it is likely to have been identified before and what you might contribute is a different method to obtain a known result.

In summary, you need the next level part of the course about translating a genetic correlation into a biological result/ mechanism. This isn't research grade work at present in my opinion.