Picard CollectGcBiasMetrics ignoring certain chromosomes/sequences What's the easiest way to run Picard GCBias ignoring certain chromosomes/sequences in the reference? Looking at the CollectGcBiasMetrics, it seems there isn't a bed file option that can be passed.

java -jar ~/picard-tools-1.134/picard.jar CollectGcBiasMetrics --help
USAGE: CollectGcBiasMetrics [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#CollectGcBiasMetrics

Tool to collect information about GC bias in the reads in a given BAM file. Computes the number of windows (of size specified by WINDOW_SIZE) in the genome at each GC% and counts the number of read starts in each GC bin.  What is output and plotted is the "normalized coverage" in each bin - i.e. the number of reads per window normalized to the average number of reads per window across the whole genome..

Version: 1.134(a7a08c474e4d99346eec7a9956a8fe71943b5d80_1434033355)

Options:

--help
-h                            Displays options specific to this tool.

--stdhelp
-H                            Displays options specific to this tool AND options common to all Picard command line
                              tools.

--version                     Displays program version.

CHART_OUTPUT=File
CHART=File                    The PDF file to render the chart to.  Required.

SUMMARY_OUTPUT=File
S=File                        The text file to write summary metrics to.  Default value: null.

WINDOW_SIZE=Integer           The size of windows on the genome that are used to bin reads.  Default value: 100. This
                              option can be set to 'null' to clear the default value.

MINIMUM_GENOME_FRACTION=DoubleFor summary metrics, exclude GC windows that include less than this fraction of the
                              genome.  Default value: 1.0E-5. This option can be set to 'null' to clear the default
                              value.

IS_BISULFITE_SEQUENCED=Boolean
BS=Boolean                    Whether the SAM or BAM file consists of bisulfite sequenced reads.  Default value: false.
                              This option can be set to 'null' to clear the default value. Possible values: {true,
                              false}

METRIC_ACCUMULATION_LEVEL=MetricAccumulationLevel
LEVEL=MetricAccumulationLevel The level(s) at which to accumulate metrics.  Default value: [ALL_READS]. This option can
                              be set to 'null' to clear the default value. Possible values: {ALL_READS, SAMPLE,
                              LIBRARY, READ_GROUP} This option may be specified 0 or more times. This option can be set
                              to 'null' to clear the default list.

INPUT=File
I=File                        Input SAM or BAM file.  Required.

OUTPUT=File
O=File                        File to write the output to.  Required.

ASSUME_SORTED=Boolean
AS=Boolean                    If true (default), then the sort order in the header file will be ignored.  Default
                              value: true. This option can be set to 'null' to clear the default value. Possible
                              values: {true, false}

STOP_AFTER=Long               Stop after processing N reads, mainly for debugging.  Default value: 0. This option can
                              be set to 'null' to clear the default value.

Given this, and wanting to calculate GCBias the using Picard (not other software), what's the best option?

The proper solution is to use a different tool (some of this you could do with computeGCBias in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s):

$ cat foo.awk
BEGIN{split(excludes, excludeList, " ")}
{exclude=0
if($1 == "@SQ") {
    for(ex in excludeList) {
        if(ex == substr($2, 4)) {
            exclude=1
            break
        }
    }
} else { \
    for(ex in excludeList) {
        if(ex == $3) {
            exclude=1
            break
        }
    }
}
if(!exclude) print}

$ excludes="chr1 chr2"
$ samtools view -h alignments.bam | awk -v excludes="$excludes" -f foo.awk | samtools -bo alignments.filtered.bam

The awk script is stripping unwanted chromosomes from the header and the body. You can then run the result through picard.


Find all the bases for given reference position Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion.
I adjusted it to find the specified position. in this case 24793.
import pysam

samfile = pysam.AlignmentFile("generated_alignment_sorted.bam", "rb" )
for pileupcolumn in samfile.pileup("chr05_modified.copy0", 10, 52000000):
    if pileupcolumn.pos == 24793:
        print ("\ncoverage at base %s = %s" %
               (pileupcolumn.pos, pileupcolumn.n))
        for pileupread in pileupcolumn.pileups:
                print ('\tbase in read %s = %s' %
                      (pileupread.alignment.query_name.split(';')[0],
                       pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

This outputs:
coverage at base 24793 = 6
        base in read m419941/6207/CCS Read=419941 = C

For all the position I tried it doesnt print out the same number of bases as it says that the coverage is.
shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?
Why is this?

pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting min_base_quality parameter when you call pileup.
To have a similar filter applied to both coverage and reads, you can use piliecolumn.get_num_aligned() to get the number of segments and pileupcolumn.pileups to get the segments themselves.


how can i download rna-seq data sample for starting rna-seq analysis by R? I am looking for an induction into RNA-seq analysis using R. I have yet to decide a data set, or a  subject area. Any help in this regard would be welcome.

You must have found the airway dataset by googling "RNASeq R tutorial"

https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html

Why won't this work for you?


Help Adjusting Bedtools Output in R I have been asked to write a script in R to edit, using data.frame, this .fa file that is the result from extracting gene sequences from a repeat masker coordinate file on the genome I am working with. Each input looks as such:

>rnd-5_family-5445_Unspecified(-)
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_Unspecified(-)
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATA
AAAATCTTTCTACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTG
TAATTTACAAAATGAATTACTTGAAGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCT
ACCGATGCCAAAAATTGACGATATTCTATTATTTAATCTAAATGTTAGTCTTTACATACCCTCCCCTA
ATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGCTAAATCTAATTCAATAGATT
GATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATTGCATATAA
CTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATG
AAAGTGATC
>rnd-5_family-2614_Unspecified(-)
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATT
AGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCT
ACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAA
ATGAATTACTT
>rnd-4_family-798_Unspecified(-)
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAA
ATAAACCTTATACTGCTAT

This is the repeat masker coordinates that the above were derived from:

>KB824701.1 417 478 rnd-5_family-5445_Unspecified . -
>KB824701.1 587 1072 rnd-5_family-2614_Unspecified . -
>KB824701.1 914 1129 rnd-5_family-2614_Unspecified . -
>KB824701.1 1138 1225 rnd-4_family-798_Unspecified . -

and ideally I would like it to be in the format of this, for example using the first one:

>rnd-5_family-5445_KB824701.1_417_478

I was told to use data.frame in R but I am not sure exactly how this will work? Can someone help please

In case an answer would still be relevant after a year and a half, below is a solution with the Biostrings package.

library(Biostrings)

# loading data.table for the fread() function.
# fread() is not only fast but is quite good at guessing how
# the file is formatted
library(data.table)

# readDNAStringSet() function creates a "DNAStringSet" object
my_fasta <- readDNAStringSet("fasta.fa")

> my_fasta
  A DNAStringSet instance of length 4
    width seq                                                                names
[1]    61 TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG      rnd-5_family-5445...
[2]   485 AATCTGAGATTAGATTTATTTATCTGTTATCC...GGCCAATCATTCAAAAGTTATGAAAGTGATC rnd-5_family-2614...
[3]   215 TTTCTAGTTTGTATTTGAGGAAAGAAAATCCT...TAGAAATTGTAATTTACAAAATGAATTACTT rnd-5_family-2614...
[4]    87 ATGCTTAATTTGCTAAGCATTTAATAACCTTT...GTAGACATTAAAATAAACCTTATACTGCTAT rnd-4_family-798_...

# the fasta header is stored at the @NAMES slot of the @ranges slot of the
# DNAStringSet object created above
> my_fasta@ranges@NAMES
[1] "rnd-5_family-5445_Unspecified(-)" "rnd-5_family-2614_Unspecified(-)"
[3] "rnd-5_family-2614_Unspecified(-)" "rnd-4_family-798_Unspecified(-)"

# removing the trailing "_Unpecified..."
my_fasta@ranges@NAMES <- sub("_Unspecified.*",
                              "",
                              my_fasta@ranges@NAMES)

my_coord <- fread("repeat_masker.txt")
> my_coord
           V1   V2   V3                            V4 V5 V6
1: KB824701.1  417  478 rnd-5_family-5445_Unspecified  .  -
2: KB824701.1  587 1072 rnd-5_family-2614_Unspecified  .  -
3: KB824701.1  914 1129 rnd-5_family-2614_Unspecified  .  -
4: KB824701.1 1138 1225  rnd-4_family-798_Unspecified  .  -
#removing the ">" sign

my_coord$V1 <- sub(">",
                   "",
                   my_coord$V1)

# appending first 3 columns of the coordinate file to "NAMES" as requested by the OP
my_fasta@ranges@NAMES <- paste(my_fasta@ranges@NAMES,
                               my_coord$V1,
                               my_coord$V2,
                               my_coord$V3,
                               sep = "_")

# writing the modified DNAStringSet object into a fasta file
writeXStringSet(my_fasta, "edited_fasta.fa")

>rnd-5_family-5445_KB824701.1_417_478
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_KB824701.1_587_1072
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTA
CCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTTGA
AGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCTACCGATGCCAAAAATTGACGATATTCTATTATTTAA
TCTAAATGTTAGTCTTTACATACCCTCCCCTAATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGC
TAAATCTAATTCAATAGATTGATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATT
GCATATAACTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATGAAAG
TGATC
>rnd-5_family-2614_KB824701.1_914_1129
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATTAGATTTATTTAT
CTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTACCACTTTTGATTCATTTTTTAAA
TCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTT
>rnd-4_family-798_KB824701.1_1138_1225
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAAATAAACCTTATA
CTGCTAT


Since they share the ensembl_ID column, you can merge them, then assign the symbol column to the rownames, then delete the symbol column.

something like:

merged <- merge(mat1, mat2, by = 0)
rownames(merge) <- merged[,'symbol']
merged[,symbol] <- NULL

https://stat.ethz.ch/R-manual/R-devel/library/base/html/merge.html


Trying to create a .bam file without the need for a .sam file I'm trying to use the code specified in this link to create a .bam file without the need for a .sam file.
Here is the code I'm using:
bwa mem -R "@RG\tID:$SAMPLE\tPL:illumina\tSM:$SAMPLE" \
    $REF_PATH/hg38.fasta \
    mapping/R1.fastq.gz \
    mapping/R2.fastq.gz | \
    samtools sort - | \
    tee mapping/$SAMPLE\.bam | \
    samtools index - mapping/$SAMPLE\.bam.bai \
    2>> log

and I get the following error:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
`[E::main_mem] fail to open file ` '.

The ref path and file locations work when creating and converting a .sam to a .bam file.
Any ideas?

This is an error from bwa saying that it can't find your input files:
Here's a checklist:

Check that you have no whitespace after your \ newline breaks. I.e. \ is the very last character
Check that you have no whitespace in any path names
Quote your variables: $REF_PATH -> "$REF_PATH"
Use set -x at the top of your bash script to debug.
In general, read this


Important genes beyond PAM50 for breast cancer classification I am currently trying to complete a Breast Cancer Classification task using Neural Networks. I have experimented with using my full dataset of gene expressions(1000+) and managed to improve my results by using only the PAM50 genes. I was wondering though what other genes I could add to the list to even better improve my results.
I have tried looking for papers about extending the PAM50, but could not find anything. I would greatly appreciate any advice.
I have also asked the question here

There is a CSV table in this paper with 33 sets of genes.


How identifiable are human omics data and how to mitigate their identifying features? Say a database were to store human omics datasets. The human subjects are known and the sample size is rather small in size initially (n=500). The database contains genomics, transcriptomics, proteomics, gut metabolomics, and epigenomics. Since the sample size is rather small initially, it would be important that certain identifying features are mitigated so that individual subjects cannot be identified. This leads me to two questions:

What identifiable features could be present in these different types of raw data? (Such as race, sex, age, hair color, eye color, height, where someone has lived, and anything else I may not even be considering)? Which of these omics types would be the most and least dangerous for identifying individuals? How can the identifiable nature of these data be mitigated?

If environmental metagenomics is also collected by these human subjects, is it possible to identify the human subjects by contamination? (i.e. some of the reads from the metagenomics data inadvertently contain human reads?) How can the identifiable nature of these data be mitigated?

I think this subject may be a bit futuristic, but I am very uninformed. If there are any references that provide additional thinking about these topics, please kindly share. Thank you for sharing your thoughts.

This is one of the major problems with genomic information in todays research. This was highlighted some years ago with the police using publicly available genomic data bases to idenitfy unkown murder suspects.
The full extent of this issue was exemplified in a
paper by Yaniv Erlich from 2018 that should be a good starting point for you.
They claim that giving the current (2018) amount genetic information more than half the searches for a person of european descent will result in a thrid cousin or better match. And these odds will only inprove with time.
Therefore, the genetic information not only can give you a hint at inherited variables like the approximate height of a person but can even yield their name and address.
About your second question, I would think that the amount of host (i.e human) material collected alongside the microbiome sample is rather small. This probably is not sufficient to obtain enough coverage for a profile. But I could be mistaken.


How do I prevent the FeatureHeatmap function from the Seurat package, from sorting my data groups in alphabetical order when plotting data? How can I prevent a function from sorting my data groups (factors) in alphabetical order without affecting the integrity of the data? I am analysing single cell RNA sequencing data using Seurat 2.3.4. I want to be able to visualise the expression of two genes, Gene1 and Gene2 across cell clusters using the FeatureHeatmap() function from Seurat. I have assigned to the metadata slot of my Seurat object the factors or conditions being tested in the experiment. However on the resulting FeatureHeatmap plot, the data are sorted and grouped in alphabetical order to my disfavour. I have attempted a solution given on stackoverflow such as the one below:

conditions <- c(rep("WTNaClctrl", 1229),
               rep("WTBlu", 1328),
               rep("ILKOBlu", 349))

conditions <- as.data.frame(conditions)

seurat_object@meta.data$conditions <- factor(conditions$conditions, levels=unique(conditions$conditions))

When I use the FeatureHeatmap() function, the plot still sorts the data groups in alphabetical order:

FeatureHeatmap(seurat_object, features.plot = c("Gene1", "Gene2"),
               pt.size = 0.25,
               group.by = "conditions",
               max.exp = 3)

This plot gets produced:

I would like groupings not to appear in alphabetical order as such:

As previously suggested StackExchange, I have tried to upgrade my Seurat object to Seurat V3:

UpdateSeuratObject(Seurat_object):
Object representation is consistent with the most current Seurat version.

devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0') : package ‘curl’ successfully unpacked and MD5 sums checked:

Error: (converted from warning) cannot remove prior installation of package ‘curl’.

Thank you in advance for any suggestions.

I had this same issue. I reordered the levels and that worked for most functions, but would switch back with the FeatureHeatmap. Turns out there is a line in the function that changes it to a character.

I ended up making a "new" function and just changed that one line from:

data.plot$ident <- as.character(x = object@ident)

to:

data.plot$ident <- object@ident

Once that was fixed, the order that I assigned worked out just fine.


Can't locate object method "new" via package Bio::Tools::Primer3 Hi there I am still a novice at this. I'm trying to run this perl script called FullSSR https://sourceforge.net/projects/fullssr, by following the instructions in their 'README.TXT' https://sourceforge.net/projects/fullssr/files/.
I did everything as instructed by the app developer(s) which included among others, replacing the Primer3.pm file provided in the *.tar.gz file into /usr/share/perl5/Bio/Tools/Run/; opening the config.txt and specifying the path of primer3_core (file) and primer3_config/ (folder); placing the sequence file I wanted to analyse into the FullSSR folder and executing the perl script using the command line: perl 'FullSSR\(ver-1.1\).p; -i sequence.fasta -m 5.
However, it was to no avail. Instead, I got the message "Can't locate object method "new" via package Bio::Tools::Primer3" (perhaps you forgot to load "Bio::Tools::Primer3?) at /usr/share/perl5/Bio/Tools/Run/Primer3.pm line 497.
I am going to try to install/load Bio::Tools::Primer3 with the cpan module (cpanm) and will see if that does it.
Perhaps someone out there could suggest a helpful (easy to understand) manual on CPAN?!

I managed to figure out the answer to my question. I installed Bio::Tools::Primer3 by typing the following into terminal: 'sudo cpanm Bio::Tools::Primer3 --force' to install it. It only works when I used 'sudo' and 'force'


What is the current state-of-the-art in assembling hybrid transcriptomes? We are considering attempting de novo assembly of a species transcriptomes (i.e. without a reference genome) using the combined NGS outputs of Iso-seq and Illumina.

One example I saw (Li et al 2017), used the standard PacBio tools to assemble a transcriptome, followed by correcting the assembled sequence using a tool called proovread, followed by cd-hit-est and then cogent.

Would this be considered a reasonable pipeline or is there an alterative recommendation?

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

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


DESeq2 complicated design - effect of replicated samples I have RNAseq data from a relatively complicated experimental design with variables = genotype, treatment, time, and batch. I have 2 biological replicates for each genotype/condition, however unfortunately in the first iteration of the experiment one sample had poor RNA quality and is unusable. To get around this lack of a replicate, we repeated the experiment for the two biological replicates in the condition of the lost sample, and did a new library prep/sequencing from this repeated experiment.  The repeated experiment obviously lead to a batch effect, and although it isn't the main driver of variance, it is a significant factor.  I would like to model or remove the batch effect from the experiment using DESeq2 if possible, but after reading this similar post: can I model technical replicates in DESeq2? I am concerned about inflating the power of the test for those conditions that were replicated.

Here is a diagram of the experimental design to make it more easy to understand:

My "colData" for the DESeq2 object looks like this:

My previous thought was to simply collapse the genotype+time+treatment variables into a single combined factor "condition" and then model the batch effect using a design of "~batch + condition".  However this doesn't really address the concern of inflating the power of the tests due to technical replication.

Should I Just collapse the counts of these replicated conditions even though the batch affect will affect the combined count values? Or should I use a different tool than DESeq2 ?  Or is there a way for me to remove the batch effect prior to collapsing the counts for replicates?

What would  you guys recommend?

I agree with @devon-ryan that the design matrix looks reasonable, i.e. that it's okay to treat the new samples as biological replicates of the original batch, and that the batch number should be included in that matrix (assuming DESeq2 doesn't complain about it) to filter out most of the batch effects.

DESeq2 tends to be conservative in its differential expression calculations. While it's true that more replicates in one sample will increase the number of differentially-expressed genes that are statistically significant (because the biological variation can be better modelled), I wouldn't expect this to change the overall picture.

However, whole-genome differential expression testing should not be used on its own to demonstrate biological relevance. I recommend as a first-pass sanity check that people compare any differential expression results to what is seen on a genome browser to identify any quirky things going on (e.g. expression across only a part of the gene). Ideally, results should be confirmed with targeted analyses on the same (or similar) biological samples, either at the transcript level (e.g. qPCR) or at the protein level (e.g. Western Blot). When aggregating evidence from multiple analyses, the specific way in which one of those analyses was done is less important.

Aside: It'd obviously be better [in the future] to make sure that your sample replicates are both robust to drop-out, and allow for outliers to be noticed, reducing the batch-effect problem. In practise this means a minimum of four replicates per informative condition, but because funders / researchers are well known for only choosing the minimum, I'd just state that six replicates per sample as what is needed to produce good results.


Changing the record id in a FASTA file using BioPython I have the following FASTA file, original.fasta:

>foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

I need to change the record id from foo to bar, so I wrote the following code:

from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        print record.id             # prints 'foo'
        if record.id == 'foo':
            record.id = 'bar'
        print record.id             # prints 'bar' as expected
        SeqIO.write(record, corrected, 'fasta')

We printed the record id before and after the change, and get the expected result. We can even doublecheck by reading in the corrected file again with BioPython and printing out the record id:

with open(corrected_file) as corrected:
    for record in SeqIO.parse(corrected, 'fasta'):
        print record.id                  # prints 'bar', as expected

However, if we open the corrected file in a text editor, we see that the record id is not bar but  bar foo:

>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

We can confirm that this is what is written to the file if we read the file using plain Python:

with open(corrected_file) as corrected:
    print corrected.readlines()[0][1:] # prints 'bar foo'

Is this a bug in BioPython? And if not, what did I do wrong and how do I change the record id in a FASTA file using BioPython?

If I use SeqIO.parse(filehandle, 'fasta') to parse a FASTA file, then it will return a SeqRecord object where the id and name are the first word (everything before the first whitespace) of the line beginning with > and the description is the complete line (all not including the initial >). (This behaviour can overruled by providing a custom title2ids function).

So, if I for example have a original.fasta file like:

>Lorem|ipsum>dolor sit amet
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

Then I will get:

>>> from Bio import SeqIO
>>> path = r'C:\path\to\original.fasta'
>>> records = SeqIO.parse(open(path), 'fasta')
>>> record = next(records)
>>> record.id
'Lorem|ipsum>dolor'
>>> record.name
'Lorem|ipsum>dolor'
>>> record.description
'Lorem|ipsum>dolor sit amet'

So, what happens in my case when you do

    if record.id == 'foo':
        record.id = 'bar'

is that the record.id is successfully changed from foo to bar, but that the record.description is not changed and stays foo. That's why, when the FASTA file is printed out, I see the described behaviour of

>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

So, the solution to my problem is to both change the id AND the description:

from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        if record.id == 'foo':
            record.id = 'bar'
            record.description = 'bar'
        SeqIO.write(record, corrected, 'fasta')

with open(corrected_file) as corrected:
    records = SeqIO.parse(corrected, 'fasta')
    for record in records:
        print record.id             # prints bar
        print record.name           # prints bar
        print record.description    # prints bar


Using preprocessing/alignment functions on the server I am new to bash and the processes behind cluster computing in general and need some help with understanding some basics.

After looking all over the internet and this forum (+ askUbuntu) I found nothing that addresses this issue.

I have a series of raw RNA sequencing files: .bam and wish to utilize samtools (or others) to begin the data preprocessing steps defined in many workflows.

I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error: -bash-samtools: command not found

When I try to install the package, say samtools, by running sudo apt install samtools, I get an error stating I do not have permission to do such.

What do I do? Must I bring up this issue to the server provider or is there a way around it?

Thank you for your time,

The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of sudo ability is not an issue.

  I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error

Installing packages on your local machine won't affect the remote server unless they share a file system, so this is expected.

As an aside, BAM files aren't really "raw", they're considered processed data (fastq/fast5 files are considered truly "raw" data in the sequencing world).


What is the fastest way to get the reverse complement of a DNA sequence in python? I am writing a python script that requires a reverse complement function to be called on DNA strings of length 1 through around length 30. Line profiling programs indicate that my functions spend a lot of time getting the reverse complements, so I am looking to optimize.

What is the fastest way to get the reverse complement of a sequence in python? I am posting my skeleton program to test different implementations below with DNA string size 17 as an example.

#!/usr/bin/env python
import random
import timeit

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

By the way, I get output like this. It varies by the call, of course!

naive implementation
       1.83880s total,
       271916.7 strings per second
       -0.7% increase over baseline
global dict implementation
       1.74645s total,
       286294.3 strings per second
       4.3% increase over baseline

Edit: Great answers, everyone! When I get a chance in a day or two I will add all of these to a test file for the final run. When I asked the question, I had not considered whether I would allow for cython or c extensions when selecting the final answer. What do you all think?

Edit 2: Here are the results of the final simulation with everyone's implementations. I am going to accept the highest scoring pure python code with no Cython/C. For my own sake I ended up using user172818's c implementation. If you feel like contributing to this in the future, check out the github page I made for this question.

the runtime of reverse complement implementations.
10000 strings and 250 repetitions
╔══════════════════════════════════════════════════════╗
║ name                       %inc   s total  str per s ║
╠══════════════════════════════════════════════════════╣
║ user172818 seqpy.c        93.7%  0.002344  4266961.4 ║
║ alexreynolds Cython (v2)  93.4%  0.002468  4051583.1 ║
║ alexreynolds Cython (v1)  90.4%  0.003596  2780512.1 ║
║ devonryan string          86.1%  0.005204  1921515.6 ║
║ jackaidley bytes          84.7%  0.005716  1749622.2 ║
║ jackaidley bytesstring    83.0%  0.006352  1574240.6 ║
║ global dict                5.4%  0.035330   283046.7 ║
║ revcomp_translateSO       45.9%  0.020202   494999.4 ║
║ string_replace            37.5%  0.023345   428364.9 ║
║ revcom from SO            28.0%  0.026904   371694.5 ║
║ naive (baseline)           1.5%  0.036804   271711.5 ║
║ lambda from SO           -39.9%  0.052246   191401.3 ║
║ biopython seq then rc    -32.0%  0.049293   202869.7 ║
╚══════════════════════════════════════════════════════╝

I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions:

import string
tab = string.maketrans("ACTG", "TGAC")

def reverse_complement_table(seq):
    return seq.translate(tab)[::-1]

The thing with hashing is that it adds a good bit of overhead for a replacement set this small.

For what it's worth, I added that to your code as "with a translation table" and here is what I got on my workstation:

global dict implementation
       1.37599s total,
       363374.8 strings per second
       3.3% increase over baseline
naive implementation
       1.44126s total,
       346919.4 strings per second
       -1.3% increase over baseline
with a translation table
       0.16780s total,
       2979755.6 strings per second
       88.2% increase over baseline

If you need python 3 rather than python 2, then substitute tab = str.maketrans("ACTG", "TGAC") for tab = string.maketrans("ACTG", "TGAC"), since maketrans is now a static method on the str type.

For those wondering, using biopython is slower for this (~50% slower than the naive implementation), presumably due to the overhead of converting the strings to Seq objects. If one were already reading sequences in using biopython, though, I wouldn't be surprised if the performance was much different.


Obtaining data table headers from GEO using GEOquery For a study in GEO, I would like to obtain the data table header descriptions, specifically the "VALUE" column for all samples in the study.

If you go here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99511)
and then scroll down and click one of the samples: let's choose "GSM2644971". Then scroll down and you should see "Data table header descriptions" and below that you should see "VALUE Normalized (provided the normalization method) Average Beta". That information is what I want.

I tried using assayData() from the Biobase package but I didn't know whether the method takes a sample, a matrix of samples, or something else as a parameter. Apparently it takes a S4 and then returns some random text but I don't know how to use that random text.

EDITED:

For example, if I do:

library(Biobase)
library(GEOquery)

getgeo<-getGEO("GSE99511")
assayData(getgeo$GSE99511_series_matrix.txt.gz)

When I use assayData() I would want it to return "Normalized (provided the normalization method) Average Beta".

Instead it returns:

<environment: 0x110674178>

Please let me know if the question doesn't make sense or if there's another method that I should be using.

The short answer is that if you are seeing information in a GEO sample, then it is the sample that you need to download and access using the appropriate method. In this case, that method is Columns():

library(GEOquery)
gsm <- getGEO("GSM2644971")
Columns(gsm)

          Column      Description
1         ID_REF
2         VALUE       Normalized (provided the normalization method) Average Beta
3         Detection   Pval

The longer answer: you will really benefit from sitting down with the GEOquery documentation and trying to understand the objects that it uses and their methods. It is not something you can work out by trial and error for every case-by-case basis. The tutorial is a really good starting point (you will find Columns() in there, for example).


gmx mdrun error: "Error in user input: Invalid input values In option s" I'm new to bioinformatics and gromacs I tried to run the code and got this error
Program:     gmx mdrun, version 2020
Source file: src/gromacs/options/options.cpp (line 179)
Function:    void gmx::internal::OptionSectionImpl::finish()

Error in user input:
Invalid input values
In option s
Required option was not provided, and the default file 'topol' does not
exist or is not accessible.
The following extensions were tried to complete the file name:
  .tpr

I already have the file topol. I searched and found that I need to run mdrun-h to ensure that the command work but when I check
mdrun -h

I got
mdrun: command not found

If the command line says "command not found", it means that the binary is not executable or not found in any of the paths mentioned in the $PATH environmental variable. You can check for example this tutorial for details on how it works and how to resolve problems related to this variable.
Which probably means, that you need to go a step back, to the installation instructions of gromacs and try to install the mdrun binary. However, if the mdrun was not installed successfully, it's very likely that more stuff of the package will be missing, so it might be worth reinstalling it all.
P.S. It's not obvious to me if getting mdrun it will resolve the first posted problem (I have never worked with gromacs). If you would share the link to the research you have done, that would help others to help you.


Link Prediction in Bipartite Networks between biological pathways and drugs Currently I'm building a big matrix (using microarray, mass spectrometry, RNAseq data) that consist by pathways (rows) and treatment/drugs (columns).

The values of that matrix are scores that describe the effect/relation of that specific treatment/drug on each pathway.

Now, some of the drugs is known to target specific disease (target disease) and thus is expected to have a high score with that disease's pathways. (A)

The rest of the drugs is known to target other diseases.(B)

My thought now is that, if I label the known drugs as "1" (or whatever), could I use a prediction model to find out new relationships (links) between drugs from the B class and the significant disease pathway from A ?

The whole thing, we can represent it as a bipartite network shown below.

For example, here the denoted as green, drugs are the known drugs belong to the class A that mentioned before and they are connected highly with some pathways of the target disease.

From the other hand blue and red ones are drugs taken from other cases (class B) that are used to target other diseases, but as my plot shows, blue seems to has an "affinity" to some of the pathways of the target disease and thus should be an alternative candidate for that disease. The red from the other side, doesn't have any "affinity" to any significant pathway and thus we can not consider it as possible candidate.

To find a solution to that kind of problem, I looked for algorithms and methods developed for social networks (e.g recommendation etc.) and mostly I read about random walks with return.

But before to proceed, I would like to know if someone else tried such an approach here and having more experience.

I am not sure if I understood correctly how you classify  the drugs. But what you attempt to do is similar to what they do here but in that article they use the drugs targeted to a molecule/pathway to find combination of drugs that are better for the disease (not an alternative drug), see the first image:

In that article they mention the pathway cross-talk inhibition (PCI), you might want to calculate it for your pathways to assess the relative reduction of the pathway efficiency.


Detecting whether genome is genetically engineered? Is there a methodology for detecting whether a genome has been tampered with by a biologist?  For instance, say there has been an ebola outbreak such as in "The Andromeda Strain".  Is there a way a researcher can tell by just looking at the ebola genetic data and determine it is a genetically engineered strain?

Some ideas I have:

organism is abnormally close to unrelated organisms, i.e. some of the genetic code from the other organisms has been inserted into the first
subsequences are very different from the rest of the genome.  If there is a consistent rate of mutation applied to the entire genome, then a subsequence that is abnormally noise free could be inserted.
combination of 1 and 2: finding clear instances of unrelated organism genomes inserted into the genome being examined.

Additionally, are there any publicly accessible reference genomes that are known to be genetically engineered, that I could use as examples to test such ideas?

UPDATE: the wikipedia article on detecting genetically modified organism is not quite what I am looking for, since it claims current methods cannot detect unknown genetic modification.

https://en.wikipedia.org/wiki/Detection_of_genetically_modified_organisms

  Currently, it is highly unlikely that the presence of unexpected or
  even unknown GMOs will be detected, since either the DNA sequence of
  the transgene or its product, the protein, must be known for
  detection.

The definition of a GMO for example UK GMO regulatory bodies (HSE) is strictly the artificial introduction of foreign genetic material using targeted genetic engineering, such as splicing,  genome integration etc... Attenuation for example by serial passage the basis of the yellow fever vaccine, does not count as a GMO.

The simple answer is look at phylogenetic trees between the GMO strain against the total distribution of isolates. If large regions have been spliced the 'silent mutation rate' looks very different,.because the GMO parent is from a very different environment to the recipient.  The switch in replication rate affects the numbers of 'silent mutation', which is fast in RNA viruses.

If the GMO is a point mutation, ie limited amino acid change then I've actually done this with dengue virus and it is fairly easy because the novel mutation make the GMO 'stick out'.of the phylogeny against all natural mutation patterns. At this point the outlier is tracked back to Genbank where the relavant publication defines the engineered virus.

If you are inferring a conspiracy theory that is a different process and a different discussion, particularly if it is Ebolavirus which is a Category 4 pathogen. The investigation is again within the phylogeny, but in this case the answer is always the isolate is a natural zoonose from a sylvatic (wild) transmission cycle which has spilled over into humans.


Knowing which positions are variants Let's say I'm sequencing a donor, and I have 20X coverage. And let's say all the reads covering a particular base has the same allele at the position. Then, is that position considered a "variant" with only 1 allele in that position or is that not even considered a variant?

The reason I'm asking is, I was given some read data, and it's in the form location:variant - location:variant.

But there was often the case where for all the variants at a particular location, there was only one allele. But this data is labelled to only have variants. Hence, how is a particular location deemed a "variant."

For example, lets say I'm told that the following is "variant reads given in the form read# location:variant - location:variant"

read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

Then, why are these two locations considered variants? It seems like it's just a heterozygous locations. Why were these locations considered to carry variants? Similarly,

read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

In those reads, it just seems like it's homozygous. Why is it being called a variant? How was it classified as such?

Put another way, if I'm given 100000 reads, and I said I only want to look at informative reads, and so I'm given a subset with only 100 reads, and those 100 reads look something like the above, then how were these reads selected?

Set 1

read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

I understand this as showing the location of 2 variants, supported by 2 reads each. Now, if the reference has an A at position 6523432, your sample could be heterozygous for an A->C variant at position 6523432. If the reference has a C at that position, then your sample could be heterozygous for C->A. If the reference has something else again, for example G, then your sample could be heterozygous for G->A/C. Similarly for location 6710432.
Set 2

read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

In this case, if the reference has an A at position 6523432, the reads do not support a variant. However, if the reference has anything else, then the reads would support a homozygous variant e.g. G->A. And the same for position 6710432.

Now, in all of these cases, you cannot know whether a variant is real without taking into account the number of reads that support it, the quality of the alignment, and the quality of the reads. That's why we do this using dedicated variant callers, software that is designed to take all these factors into account and decide whether to call a variant at that position.

But in all cases, the variant must be defined with respect to a reference genome. You cannot call a variant just by looking at reads without any knowledge of a reference. A "variant" just means "not what is found in the reference". It is largely a matter of random chance whether something is a "variant" or "normal". If the individual(s) whose genome was taken as a reference has the same sequence, then you call it normal. If not, then you call it a variant. There are many cases where the standard hg19 or hg38 genomes have a "normal" residue which is actually far less common than the known "variant". There's nothing special about the reference genome, after all, it's just the one we chose to call a reference.


What is "aligned sequences" and "consensus sequence" in the context of sequence logo? How to compute these? In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).

so, What is "aligned sequences" and "consensus sequence"?

how can I get these sequences based on following data?

Figure 2.5 from the book "Machine Learning - A Probabilistic Perspective"

more info about this figure at here

Panel (a) shows a sequence gap-free alignment. Each row corresponds to a contiguous 15 base pair sequence of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.). Each column corresponds to a particular position/residue (e.g. column 1 is the 1st position of the alignment etc.).

Panel (b) shows the sequence logo. As noted in the comments, this does not correspond to the alignment in (a). Instead, here is the corresponding logo:

From this you can intuitively understand that at the 1st sequence position you have only As and Ts and their frequency corresponds to the y-axis. For the 3rd position, you have only As, so only A is represented here in the logo etc.

The consensus sequence is the most frequent letter at each position in the alignment, i.e. [AT]TAGC etc. where [AT] means A or T since they are jointly high in frequency (5 each).


extract chimeric and multimap reads from bam file I am trying to extract all chimeric and multi-map reads from either SAM/BAM file. Is there any simple command to do that? Can I use htslib for parsing sam/bam files and extracting based on any flags. Any help on how to start will be appreciated.

This is how I do it.
samtools view -@5 -f 0x800 -hb /path/sample.bam > sup.bam

Then if you want it as a fasta
samtools view -F 0x1 -hb sup.bam | samtools fasta -F 0x1 - > sup.fa


Comparing groups of 16S rRNA sequences - how to? I have several groups of 16S rRNA sequences associated with taxonomic groups (let's call them A, B and C). The sequences within each group have a common ancestor and are on average more closely related to sequences within the same group. The groups are also related, but less.
Ideally I'd like to have a table like this (where percentages represent the average similarity):
           Group A     Group B     Group C
Group A      99%
Group B      97%         98%
Group C      94%         95%         99%

Is there software that can calculate these percentages for me? I'm thinking that I could set up a many-against-many BLAST search and then try to use pandas in python in order to parse the resulting table...But if something like this already exists, I'd rather not reinvent the wheel. Thank you!

If you want to make an all vs all with blast this can be helpful. You can create a blast database with your sequences and then search against itself.
makeblastdb -in 16S_sequences.fasta -dbtype nucl -out my_16S_sequences_db

blastn -db my_16S_sequences_db -query 16S_sequences.fasta -outfmt 6 \
       -out allvsall_16S.tsv -num_threads "$(nproc)"

If you want a percent identity matrix for a pairwise alignment you can try with this approach: https://www.biostars.org/p/220154/


Are there any databases of templates for common bioinformatic file formats? I want some templates of different file formats that I can use to test my scripts and identify possible bugs in my code.

For example, consider nucleotide FASTA, a simple but often abused format, I would want templates to capture regular and irregular formats, like I have seen all of these:

1) Single line sequence

>1
ATG

2) Multi-line sequence

>1
AT
G

3) Upper and lower case letters in sequence

>1
Atg

4) Ns and Xs (and possibly other letters) in sequence

>1
ANnxX

5) Unusual headers (sometimes non-ASCI characters, need to consider the encoding)

>ATG >汉字
ATG

6) Whitespace between records

>1
ATG

>2
ATG

7) Duplicated headers

>1
ATG
>1
ATC

8) Empty headers or sequences (valid FASTA?)

>
>

9) No new line '\n' character on last line (can mess-up file concatenation)

>1
A# < no new-line here

10) Different newline characters depending on the OS

>1
A# \r\n vs \n

etc.

There should be separate templates for nucleotide and protein FASTA, and separate ones for aligned FASTA.

It would ideally include other aspects too, like different compression formats (such as .gz, .bzip2) and different file extensions (such as .fa, .fasta).

I have never seen resources that provide templates covering these, but I think it would be useful. Of course I could build my own templates but it would take time to capture all the likely variations of the formats, particularly for more complex file formats.

Note, I am not just interested in FASTA format, it was an example.

Also note, I know about tools (such as BioPython) that should handle many formats well, but they may also have bugs. Anyway, in practice sometimes I end up parsing files myself directly because I don't want the overhead or dependency of an external package.

EDIT: Please don't answer this question to say you don't know of any such resources, me neither, hence the question. bli's helpful answer shows that there is at least one test suite that could be used as a starting point. I know that it is normally easy to look up the specification of any particular file format.

You mention Biopython, which contains tests: https://github.com/biopython/biopython/tree/master/Tests.

Some of the tests consist in reading files present in the folders listed in the above link. These files could be a starting point for a database of test files. Whenever one comes across a test case not covered with these files, one could construct a new test file and contribute it to Biopython, along with a test, or at least file an issue: https://github.com/biopython/biopython/issues.

That would be a way to contribute to Biopython while constituting a database of test files.


Changing a wide range of colours to a limited gradient I returned a FeaturePlot from Seurat to ggplot. My plot has a weird range of colours as below

I produced this plot by this code

> head(mat[1:4,1:4])
             s1.1        s1.2 s1.3       s1.4
DDB_G0267178    0 0.009263254    0 0.01286397
DDB_G0267180    0 0.000000000    0 0.00000000
DDB_G0267182    0 0.000000000    0 0.03810585
DDB_G0267184    0 0.000000000    0 0.00000000
>

I have converted expression matrix to a binary matrix by 2 as a threshold

mat[mat < 2] <- 0
mat[mat > 2] <- 1

> head(exp[1:4,1:4])
             s1.1 s1.2 s1.3 s1.4
DDB_G0267382    0    0    0    1
DDB_G0267438    0    0    0    1
DDB_G0267466    0    0    0    0
DDB_G0267476    0    0    1    0
>
> exp=colSums(exp)
> exp=as.matrix(exp)
> colnames(exp)="value"
> exp=as.data.frame(exp)
> cc <- AddMetaData(object = seurat, metadata = exp)
> cc <- SetAllIdent(object = cc, id = "value")
> TSNEPlot(object = cc, do.return= TRUE)

How I can convert this range to a gradient of colours for example in 8-18, 18-28, 28-38, 38-48 range in blue to yellow please? Something like below

Thank you for any help

Then by ggplot now I scaled my colours but I don't like my clusters as so and I don't know how to retain my clusters as a featureplot by this new color gradient

> head(cc@meta.data)
     nGene    nUMI    orig.ident res.0.7 CELL STAGE GENO dataset stage.nice celltype value
s1.1  4331  373762 SeuratProject       0 s1.1   H16   WT       1        H16        0    34
s1.2  5603 1074639 SeuratProject       0 s1.2   H16   WT       1        H16        0    26
s1.3  2064   49544 SeuratProject       0 s1.3   H16   WT       1        H16        0    27
s1.4  4680  772399 SeuratProject       1 s1.4   H16   WT       1        H16        1    29
s1.5  3876  272356 SeuratProject       1 s1.5   H16   WT       1        H16        1    21
s1.6  2557  122314 SeuratProject       0 s1.6   H16   WT       1        H16        0    31

> ggplot(as.data.frame(cc@meta.data), aes(x = cc@meta.data$CELL, y = cc@meta.data$res.0.7, colour =cc@meta.data$value)) +
+     geom_point(size = 5) +
+     scale_colour_gradient(low = "yellow", high = "blue")

By below code I obtained a tsne in link

> cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))
>
> TSNEPlot(cc, colors.use=cols)

enter link description here

Now I want to know how I could convert this range to a 8-18, 18-28, 28-38, 38-48 colour range as a gradient of blue to yellow?

TSNEPlot()

TSNEPlot() will always treat your variables as discrete. My approach is to manually generate a gradient with unique colors for each factor level and pass it to the cols.use argument in TSNEPlot().

#generate values for testing purposes, one value for each cell
value <- sample(seq(from=8, to=48, by=1), size = length(rownames(unfiltered_cca@meta.data)), replace=TRUE)
names(value) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")
exp <-  SetIdent(exp, ident.use=exp@meta.data$value)
#create a gradient between specified colors, multiply by sequence to get appropriate length
#length.out is the number of different colors, number of factor levels
cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))

TSNEPlot(exp, colors.use=cols)

I would present this with no legend:

TSNEPlot(exp, colors.use=cols, no.legend = TRUE)

FeaturePlot()

You can also simply use FeaturePlot() instead of TSNEPlot() to visualize the gradient. Using the same data as above:

FeaturePlot(object = exp, features.plot = "value", reduction.use = "tsne", no.legend = FALSE, cols.use = c("beige", "red"))

You ask for a continuous scale, but this is not what is shown in your second plot.  You have a different color for each discrete range of values. I assume here that you want a continuous gradient. To achieve something like the second plot you provided requires a different approach.


Finding the location and unit length of repetitive sequences within a long read After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from Nippostrongylus brasiliensis [each base is colour-coded as A: green; C: blue; G: yellow; T: red]:

[FASTA file associated with this sequence can be found here]

These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:

Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times
Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal

After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.

For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:

Use LAST to produce a dot plot of self-mapping for the mapping
Visually identify the repetitive region, and extract out the region from the sequence
Use a combination of fold -w <width> and less -S to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size
Display the sequence in a rectangular and circular fashion using my own script, wrapping at the repeat unit length

But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.

Is there any better way to do this? Given an arbitrary DNA sequence of length >10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?

An example sequence can be found here, which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.

A Kmer-based Analysis

I've now seen human sequences with repetitive regions in excess of 10kb (i.e. out of the range of most linked-read applications). My current idea is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:

Process the sequence, storing the location of each kmer
For each kmer, find out how many times it appears in the sequence
For repeated kmers, find out how much of a gap there is between the next time that kmer occurs
Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence

Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).

It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?


how to do analysis of a table content SNP? I have a table with information about SNPs, the table has a total of 55127 rows.

How can I know the number of gene mutations in the table and if a gene is present?

The table is in xls format.

For working with a relatively small dataset in tabular format, I would recommend dplyr. You can read in a excel file using readxl

library(readxl)
library(dplyr)

df <- read_excel(FILENAME)

To get the number of mutations per gene:

df %>% group_by(`Gene name`) %>%
  summarise(num_mutations=n())

If you want to know if there are mutations in a particular gene:

df %>% filter(`Gene name` == QUERY)


Is there a convention or standard of the highest bits in a sequence logo? In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).

wiki says

  The total height of the letters depicts the information content of the
  position, in bits.

and then gives this example

  A sequence logo showing the most conserved bases around the initiation
  codon from all human mRNAs (Kozak consensus sequence). Note that the
  initiation codon is not drawn to scale, otherwise the letters AUG
  would each have a height of 2 bits.

it seems that there is a 0.1 bit letter.

is this in the context of 2bits? is 2bits the highest bits in this sequence logo? Is there a convention or standard of the highest bits in a sequence logo?

This is explained in the original paper [1] describing the logos (emphasis mine):

  The importance of a particular position in a binding site is more clearly
  and consistently given by the information required to describe the pattern
  there [7, 8]. To choose one symbol or state from two equally likely possibilities
  requires one "bit" of information. For example, to communicate the result of
  a coin-ip to someone requires 1 bit of information because only one yes-no
  question needs to be answered: "Is it heads?". If a position in a binding
  site always contains one base (e.g. G), then we need exactly two bits of
  information since two yes-no questions need to be answered: "Is it A or G?"
  (i.e. is it a purine?) and "Is it A or C?". (If the answers to both questions
  are "no", it must be a T.) Furthermore, if a position contains two bases (e.g.
  sometimes A and sometimes G), only one question suffices since a two out of
  four choice is equivalent to a one out of two choice. Therefore, only one bit
  is needed to describe a position in a binding site that contains only purines,
  but two bits are needed to describe a position that always contains adenine.

So yes, 2bits is the highest value you will see on the y-axis of a DNA sequence logo (proteins are a different story since there are more possibilities). In the one you show, the scale is from 0 to 0.2, but it's the same principle.

However, you should think of logos as a useful visual aid, and not a way of getting rigorous mathematical information about a sequence. So don't put too much weight on the bit values. If you need that sort of data, calculate it from the alignment directly.

[1] Schneider TD, Stephens RM., Sequence logos: a new way to display consensus sequences., Nucleic Acids Res. 1990 Oct 25;18(20):6097-100.


Python package or CLI tool outputting mutations of sequence with respect to reference I'm new to bioinformatics, and I'm looking for a tool that takes two FASTA files, one containing a reference, let's call it A.fasta, the other another sequence, let's call it B.fasta, and outputs just the mutations that B contains with respect to A.
I'm open as to the output format. Is there a standard for this?
Background: I'm trying to analyse Sars-CoV-2 sequences and identify mutations with respect to the reference. I've done a google search for a tool that does what I want to but couldn't find something straightaway. Is there no standard tool for this?

If it's just a comparison of two sequences that you're looking for, then diffseq from the emboss toolkit should be okay:
https://www.bioinformatics.nl/cgi-bin/emboss/diffseq
I think this tool requires the sequences to be oriented in the same direction. If you discover a lot more variants than you expect, try reverse-complementing one of the sequences.


How to select high quality structures from the Protein Data Bank? Models of structures deposited in the Protein Data Bank vary in the quality, depending both on the data quality and expertise and patience of the person who built the model. Is there a well-accepted subset of the PDB entries that has only "high quality" structures? Ideally these structures would be representative for classes of proteins in the whole PDB.

based on a real question from biology.SE

There is a very nice database, pdbcull (also known as the PISCES server in the literature). It filters the PDB for high resolution and reduced sequence identity. It also seems to be updated regularly. Depending on the cut-offs, you get between 3000 and 35000 structures.

If you are specifically interested in rotamers, you may want to look at top8000 instead, where they have checked for high resolution, and good MolProbity scores. They also provide a rotamer database.

PDB also provides their own clustering. They first cluster the sequences, and then extract a representative structure for each one, based on the quality factor (1/resolution - R_value). This has the advantage of being comprehensive, but you will have bad structures when no good ones were ever obtained.


Running tophat dockerfile in background I need to use tophat2 to align some sequences, and I wish to use the docker container. These are large sequences so it will take a long time, plus I'm working on a university server so the chances of being kicked off the server are very high. How can I use nohup and & to run tophat2 through docker, if I need some user interaction?

Assuming your user is in the group permitted to use docker:

docker pull quay.io/biocontainers/tophat:2.1.1--py27_3
nohup docker run -v /some/path:/data quay.io/biocontainers/tophat:2.1.1--py27_3 tophat2 ... &

The only thing to remember is to mount everything you'll need into the docker container (/data in the example above). If you need to debug things then you can add the -it flag to docker run and use bash as the command.

As an aside, you might find using screen or tmux preferable to nohup ... &, since you can reconnect to sessions with them.


Converting a BAM file into VCF I have NGS illumina RNA-seq reads from M. musculus (mm10). I am trying to find variants along the strand portion of the reads in the refseq (mm10).

I mapped a pair of sequence files and generated a BAM file. Now I need to generate a VCF file from this BAM file. I’ve checked Galaxy for direct tools, but I couldn’t identify one. Any suggestions on how to turn a BAM file into VCF? Thank you in advance.

It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing.

How to call variants (a vcf file) from mapped reads (a bam file) is very broad question, it depends on the sequencing technology (Illumina, PacBio, ...), the type of sequencing (whole genome sequencing, RAD, RNA-seq, ...), the sampling (calling of individuals vs population variability), the quality of a reference (human vs denovo assembly made by a postdoc in your lab)... For a lot of cases the answer is simply GATK, partially because it has well written manual, but mainly just because many people use it.

Some more detailed cases are well answered in more specific questions :

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
Somatic tumor only variant calling?
Why are there missing calls in a VCF file from exome sequencing?


Why total RNA-seq usually yields low mapping rate? Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)?

Doesn't the genome fasta file contain ribosomal RNAs (which are abundantly expected in the total RNA-seq libraries) as well?

I don't have the numbers for human right now, but I recall previously having low mapping rates for human data as well. What triggered me at the moment is actually when I mapped zebrafish data. I got around 60% mapping to the genome and 45% to the transcriptome with one dataset and 36% - 35% with another dataset (note that these are from 2 different studies and both are total RNA-seq, mapping was done using STAR and Salmon, which have their own methods).

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (this behavior can be changed with the --outFilterMultimapNmax option).
To confirm if this is actually the case, you could check the number of multi-mapping reads in the aligner's log files.


Raw Data frame manipulation in python Using python 3 I need to process qPCR sequencing raw data outputs by searching for the first occurrence of a user defined string and then making a new data frame using all lines after that string. I am trying to find solutions in the pandas doc but so far unsuccessful.

This is a raw output .csv file that I need to process. (couldn't paste complete csv as exceeds character limit, this is lines 40-50 and am hoping this text is useful?). I need to tell pandas to create a new data frame that 1. starts at the line containg the first occurance of str("Sample Name") with that line as header and containing all lines following. And then 2., only including columns ("Sample Name"), ("Target Name"), ("CT").

Could someone please help me so that I can use python to help me analyze biological data?
Many thanks,
Luke

40,Quantification Cycle Method,Ct,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

41,Signal Smoothing On,true,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

42,Stage where Melt Analysis is performed,Stage3,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

43,Stage/ Cycle where Ct Analysis is performed,"Stage2, Step2",,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

44,User Name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

45,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

46,Well,Well Position,Omit,Sample Name,Target Name,Task,Reporter,Quencher,Quantity,Quantity Mean,SE,RQ,RQ Min,RQ Max,CT,Ct Mean,Ct SD,Delta Ct,Delta Ct Mean,Delta Ct SD,Delta Ct SE,Delta Delta Ct,Automatic Ct Threshold,Ct Threshold,Automatic Baseline,Baseline Start,Baseline End,Amp Status,Comments,Cq Conf,CQCONF,HIGHSD,OUTLIERRG,Tm1,Tm2,Tm3,Tm4

47,1,A1,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.357698440551758,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9588544573203085,N,Y,N,81.40960693359375,,,

48,2,A2,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,24.05980110168457,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,15,Amp,,0.9592687354496955,N,Y,N,81.40960693359375,,,

49,3,A3,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.012556076049805,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,16,Amp,,0.9592714462250367,N,Y,N,81.40960693359375,,,

50,4,A4,False,fla11fla12-1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.803699493408203,24.419523239135742,0.5669151544570923,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9671570584141241,N,Y,N,81.40960693359375,,,

This is the code that I have so far:

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

This doesn't make much sense ... you read in Excel to a dataframe, you then read out to csv - the dataframe doesn't have an index, but you include an index and then read it in again. This could do more harm than good.

Anyway your supervisor is right it is very easy to parse this data using pandas, although there will be more expertise doing this via R.

I don't want to peer through the gloom of a CSV output, but I would like to highlight how easy it is to parse this stuff into a dataframe. The dataframe 'data' is fine.

Parse You can simply remove the top rows of bumf data = data.iloc[6:] (I might have miss counted, plus see below),
Headers You will probably have to remove the "headers" row to a list and
then introduce back onto the dataframe using data.columns =
['header1', 'header2' ...] ... rather this is the list (array in Perl) which you can automatically generate using a magic to_list type pandas command which I forget, its not list - its cuter
Columns You next might use drop to remove the columns you don't want (that might be alot of columns) or else make a subset of
the dataframe into a new dataframe, you can also use iloc for
columns too.
Index You then assign A1, A2, ... etc  as an index data.index = ['A1',
'A2', ]...etc..  you can reshape column 3 to automate this
to_numeric You might have to assign the CT values as floats
    or better to_numeric
Extrapolation You would then define a subroutine (or whatever they call them in
Python) so perform the extrapolation and use the apply command on
your dataframe and assign it to a new column, where the input is your CT values.
Output Finally, issue your to_csv command to_pickle is also a good idea because you have a df that is worth importing at a future date.
Plotting this is easy using data.plot which appears to be some sort of inbuilt matplotlib, keep in mind seaborn DONT try using R's ggplot2 within Python using rpy2 .... its a NIGHTMARE to get working and doesn't support the latest ggplot2 (!!!). R
s ggplot2 is stronger  than seaborn or matplotlib (providing you are working within R), but these two are good for quick graphs.

If I had time I'd write the code for you, but I'm afraid I'll have to leave it there. I'm still in debt with a Perl efetch script I promised a previous OP.

In summary, Python pandas is fabulous stuff for this purpose, I love it. Doing this via Perl and "subprocessing" inside R... its doable, well in this particular case I suspect its just as easy, but generally its alot more code and often fiddly.


Edit all the fasta headers using awk kindly your comment highly appreciated

have more than 4k header sequences look like:

>LTR_retrotransposon100_Gypsy?
cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

>LTR_retrotransposon10011_ERVK
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon10011_ERVL-MaLR
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon100_Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

I would like to use awk or sed to remove everything between the underscore, including the underscore and replace it with /
plus add sequential RNA followed by #.

and the target output could be:

>RNA1#LTR/Gypsy?
ctcagcagcactatttacaatagccaagacatggaaacaacctaaatgtcttatcaatag atgactggataaaggagctgtggtatatctatacaatggaataccatcagccataaaaaa gaataaaatattgccatttgcagcaacatggatggacctggagattatcattctaaggga agtaagccagaaagagaaagaaaaataccatatgatatcacttatatgtggaggtaaaaa aaaaaaaaaagacacaaatcaatttatttgcaaaacatacaTGGActttcagacatagaa

>RNA2#LTR/ERVL-MaLR
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA3#LTR/ERVK
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA4#LTR/Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

EDIT after OP's update:

Try this Perl one-liner:

perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file

Here, the command line flags are:
-e: tells the perl interpreter to look for the code inline, rather than in a file with the script,
-p: loop over the input one line at a time, execute the code in the one-liner for each iteration of the loop, and print the result into STDOUT.
The code has regex replacement (see perlre for more details):
s{>([^_]+)_.*_}{>RNA${i}#${1}/}:
replace an underscore, followed by zero or more occurrences (* quantifier) of any character (.), followed by an underscore, with slash (/).
Keep > and the original header part before the underscore, captured by parens into $1 variable.
Add string RNA${i}, where $i is the counter of the fasta headers. The counter is incremented by 1 when the substitution occurs, that is, when the next header occurs.
The variables are used like so in the substitution: ${i}, ${1} to prevent wrong interpolation.
I am also using s{}{} as delimiters for readability (as in most cases) instead of the default s///.

EXAMPLE:

echo ">LTR_retrotransposon100_Gypsy?\nACGT\n>LTR_retrotransposon10011_ERVK\nTCGA" | perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}'

Prints:

>RNA1#LTR/Gypsy?
ACGT
>RNA2#LTR/ERVK
TCGA


Compare clusters from different datasets I have two scRNAseq datasets, and I select a cluster from one of them, and another cluster from the other. I want to find differentially expressed genes between the two clusters. What would be the right way of doing that? Right now I am using DESeq2 in the following way. I select only the selected clusters' cells, merge them together forming a matrix and then feed them into DESeq2:

mtx <- merge(data_1,data_2,by="row.names",all=TRUE)
mtx[is.na(mtx)] <- 0
mtx_conv <- matrix(as.numeric(unlist(mtx[,2:ncol(mtx)])),nrow=nrow(mtx))
mtx_int <- round(mtx_conv, 0)
rownames(mtx_int) <- as.character(mtx$Row.names)

cluster <- as.factor(c(rep(0, dim(data_1)[2]), rep(1, dim(data_2)[2])))
colData <- data.frame(
  Cluster = cluster
)

dds <- DESeqDataSetFromMatrix(countData = mtx_int,
                          colData = colData,
                          design = ~ Cluster)
dds <- DESeq(dds)

I am not sure though that it is the correct way to go because I am not using the full original datasets. The reason is that I want the computations to run faster. Another way of doing that is probably combining the whole two datasets, defining 3 clusters: two groups I am interested in, and all other cells, and then running the analysis. Would it differ significantly? Am I on the right track at all?

I have another idea: to use Seurat package. There, I would run CCA algorithm to align two full datasets, and then run FindMarkers function between the two clusters. I feel like it may be wrong, because the two datasets may need to be re-normalized together but Seurat does not seem to be doing that:

https://satijalab.org/seurat/Seurat_AlignmentTutorial.html

Any suggestions would be greatly appreciated.

I've been working with the multi-CCA alignment workflow in the Seurat package and am very impressed with functionality. As usual, the CCA tutorial is excellent for an alignment workflow and will get you off the ground quickly. You can even use the DESeq2 algorithm for determining differential expression from within Seurat.

When you use the function RunCCA() you can pass rescale.groups=TRUE if you wish to re-scale.

Abstracted workflow

Data preprocessing and gene selection
Define a shared correlation space with canonical correlation analysis
Remove dataset specific cells
Align correlated subspaces using dynamic time warping
Integrated analysis across datasets (clustering, trajectory building, etc ...)

Butler, A., and Satija, R. (2017). Integrated analysis of single cell transcriptomic data across conditions, technologies, and species.


gprofiler output ENSG to gene name I would like to ask if there is any tool or R script that translates Ensembl IDs to gene names in a gprofiler GO output. An example of this file is the following one:

This is usually done using AnnotationDbi in combination with the annotation database for human org.Hs.eg.db.
library(org.Hs.eg.db)
genes <- c("ENSG00000102055", "ENSG00000232605", "ENSG00000136828",
           "ENSG00000159197", "ENSG00000201394", "ENSG00000278580",
           "ENSG00000213440", "ENSG00000166816", "ENSG00000275895",
           "ENSG00000105723")
select(org.Hs.eg.db, genes, c("ENSEMBL", "ALIAS"), "ENSEMBL")

This returns:
           ENSEMBL     ALIAS
1  ENSG00000102055       I-4
2  ENSG00000102055  PPP1R2P9
3  ENSG00000102055   PPP1R2C
4  ENSG00000232605      <NA>
5  ENSG00000136828   RALGEF2
6  ENSG00000136828  RALGPS1A
7  ENSG00000136828   RALGPS1
8  ENSG00000159197     ATFB4
9  ENSG00000159197      LQT5
10 ENSG00000159197      LQT6
11 ENSG00000159197     MIRP1
12 ENSG00000159197     KCNE2
13 ENSG00000201394   RN5S140
14 ENSG00000201394 RNA5SP140
15 ENSG00000278580      <NA>
16 ENSG00000213440      <NA>
17 ENSG00000166816     DLACD
18 ENSG00000166816       DLD
19 ENSG00000166816      LDHD
20 ENSG00000275895   U2AF1L5
21 ENSG00000105723     GSK3A

If you want to apply a 1:1 mapping of gene IDs to common names, you can use the "SYMBOL" column. In a data.frame, you can do something like this:
ensembl_to_symbol <- function(e) {
  hits <- select(org.Hs.eg.db, e, "SYMBOL", "ENSEMBL")$SYMBOL
  # Gene IDs without symbols should just stay as IDs
  no_hit <- is.na(hits)
  hits[no_hit] <- e[no_hit]
  hits
}

recode_symbols <- function(s) {
  e <- ensembl_to_symbol(unlist(strsplit(s, split = ",")))
  paste(e, collapse = ",")
}

df <- data.frame(
  genes = c(paste(genes[1:5], collapse = ","),
            paste(genes[6:10], collapse = ",")
  )
)

df$genes <- sapply(df$genes, recode_symbols)

This isn't the most vectorized (fast) solution, but it's probably easier to understand this way.


How to calculate the number of peaks that are upstream/downstream of some other peaks I have 3 histone marks,I have used Macs2 for peak calling and diffBind to analyze the peaks. I was wondering if you know any way to calculate the peak numbers of one specific histone mark that occur for example within 2kb of peaks belong to another histone mark.
Much appreciated

You can use windowBed from bedtools with -c if you just want counts. Or use -wa -wb and then use groupBy for any other operations.


How can I find the chromosomal location of a list of genes? I have a list of genes nearly 20000:
gene name (column1) and coordinates (columns 2 and 3)

head genes
    Ma0_73Ns    2359    2923
    Ma0_73Ns    4727    4868

I want to find the chromosomal location of each gene, on which chromosome each gene located? my study organism is sunflower

Ok so I found two websites that could help. The one I found first is more intuitive to go through, but it looks like it could be outdated as the second website I found seems more recent so it probably has a more complete genome to use as the reference. (Sorry, I haven't worked personally with sunflowers and don't know for sure where you would go for the "official" reference genome, or if there even is an "official" one, but hopefully one of these websites will help.)

1) This is a website I found that looks like it is a website for the sunflower genome and includes gene annotation information. (The recent publications list only goes up till 2015 though so it may not be the most recent data)

https://www.sunflowergenome.org/annotations-data/

If you scroll down to the bottom of this page there is a text file detailing the descriptions of how the files are formatted.

For example, the linked to file gives the following:

  Ha412v1r1_CDS_v1.0.fasta.gz

  Description: Nucleotide CDS for each gene
  Header format: >TranscriptID GeneID Chromosome_start_end_strand
  Example: >Ha1_00044280-RA Ha1_00044280 Ha1_94590914-94591860_-

If you download this fasta file you could then create a python script that will read in the file and extract the headers. If each header is a string variable called "mytext", you could then extract the TranscriptID, GeneID, and chromosomal info and using mytext.split(" "). You could then split on "_" to further parse the chromosomal location info.  You could save the extracted information in a dictionary using the gene name as a key and parse through the original file you mentioned containing the list of genes, in order to get the chromosomal location information you desire.  The one caveat is that this is assuming that your gene names are the same as the geneIDs the above website uses.

I did not exhaustively examine the list of files and perhaps there is another one that would be a better one to parse through for your purpose, but the linked to text file gives more details on the files and their formats and presumably they could be parsed in a similar manner.

2) Double checking my first answer I found this article:

  Badouin H, Gouzy J, Grassa CJ, Murat F, Staton SE, Cottret L,
  Lelandais-Brière C, Owens GL, Carrère S, Mayjonade B, Legrand L. The
  sunflower genome provides insights into oil metabolism, flowering and
  Asterid evolution. Nature. 2017 Jun;546(7656):148.

Which links to a website as follows:

  "The integrative web interface Heliagene provides visualization, querying tools for data mining and network exploration for the community (https://www.heliagene.org)."

So it looks like this website also has reference genomes for sunflowers, although a lot of the annotation information seems to be locked unless you have a login. Like I mentioned previously, this website is a little less intuitive to go through but you might be able to find something helpful here if you think the other website isn't official or recent enough.


Subset data frame by similar elements from rows and column between two data frames My first data frame whose first column i intend to use which is PATIENT_ID

PATIENT_ID  SEX RACE    FAB AGE BM_BLAST_PERCENTAGE WBC
TCGA-AB-2805    Male    WHITE   M0  77  67  92
TCGA-AB-2822    Male    WHITE   M0  65  99  2.9
TCGA-AB-2831    Male    WHITE   M0  59  61  1.4
TCGA-AB-2890    Male    NH/C    M0  69  60  6.1
TCGA-AB-2907    Male    WHITE   M0  69  82  1.6
TCGA-AB-2913    Male    WHITE   M0  61  95  61.6
TCGA-AB-2922    Male    WHITE   M0  83  75  77.9
TCGA-AB-2936    Female  NH/C    M0  60  50  2.1

My second data frame where i have the respective patient ID as individual column with expression values

Hugo_Symbol TCGA-AB-2803    TCGA-AB-2805    TCGA-AB-2806    TCGA-AB-2807    TCGA-AB-2808
A1BG    0   0   0   0   0
A1CF    0   0   0   0   0
A2BP1   71.4095 48.5475 50.8339 45.7498 47.1328
A2LD1   2.0576  3.3467  5.7288  0.993   5.4988
A2ML1   467.0782    16.0643 7.0019  343.5948    255.9152
A2M 11.3169 0.6693  17.823  116.1867    182.2467
A4GALT  0   2.6774  1.2731  4.9652  0.7855
A4GNT   0   0   0   0   0
AAA1    514.4033    532.7979    504.1375    697.1202    443.8335
AAAS    0   0   0   0   0
AACSL   190.3292    74.2972 84.6595 559.0864    317.3606
AACS    0   0   0   1.9861  0
AADACL2 0   0   0   0   0
AADACL3 0   0   0   0   0
AADACL4 1.0288  0   0.6365  0   0

What I tried so far transposing the second data frame so that i can get the individual column as my first column ,and subset by using dplyr library but so far im not able to ..any help or suggestion how to do it would be highly appreciable

I would like to have my expected data frame as such

PATIENT_ID  Gene B  GeneA   GeneC   GeneA   GeneA
TCGA-AB-2803-03 0   0   0   92.5926 876.5432
TCGA-AB-2805-03 0.6693  1.3387  0   92.6305 421.6867
TCGA-AB-2806-03 0   14.1884 16.3654 119.599 625.0796
TCGA-AB-2807-03 2.006   9.3247  13.5154 178.4211    537.2393
TCGA-AB-2808-03 0   42.3016 28.3975 170.2042    653.5742
TCGA-AB-2810-03 0   1.5668  8.9665  152.0747    1087.5648
TCGA-AB-2811-03 0   0   6.9556  30.3103 1414.6602
TCGA-AB-2812-03 0   14.3289 11.9466 120.3591    962.5575
TCGA-AB-2813-03 0   3.4938  1.3813  51.1883 594.7593
TCGA-AB-2814-03 0   56.8659 118.7051    87.1817 472.6912
TCGA-AB-2815-03 0   2.0692  3.6986  104.6575    502.5234

You can take the PATIENT_ID from the first dataframe and use them to select from your second dataframe. Than transpose, so something like this:

ID_Names <- First_df$PATIENT_ID

rownames(Second_df) <- Second_df$Hugo_Symbol

Selected_df <- Second_df[, ID_Names]

Expected_df <- t(Selected_df)


Converting Ensembl Gene IDs to Entrez Gene IDs through biomart Well, I'm trying to convert a list of Human Gene referenced by Ensembl Gene IDs to Entrez Gene IDs. I have been advised to use biomart.

I tried to get a kind of conversion table for all human genes. I don't know if my settings are wrong, but I didn't find any checkbox for Entrez Gene Id in the section Attributes > External References. I've just found the Entrez transcript name ID (see below), but it's not what I need.

How can I use biomart for such a conversion?

Edit:

According to this video tutorial EntrezGene ID should be an option in biomart (see screenshot below).

That video tutorial is outdated.

The ID you need is the NCBI gene ID, which is the same as the EntrezGene ID.


Fitting a principal curve into a diffusion map I returned a Seurat object with a diffusion map by

 seurat_object <- RunDiffusion(seurat_object,genes.use = seurat_object@var.genes)

as below

Now, how I could know the order of a cell's projection onto this principal curve representing its predicted progression through the development by fitting a principal curve on my diffusion map?

I know Seurat slot for DM1 and DM2 is inside the seurat@dr

by

fit=principal_curve(as.matrix(seurat_object@dr$dm@cell.embeddings))
plot(fit)
line(fit)
points(fit)

I draw this graph of my cells

But how I know the first and last cells on curve? I mean how I know cells from each cluster is being placed where on the curve? By fit$ord I will obtain orders of cells on the line but how I know the first and end of the line?????

If fit$ord gives you the order of cells on the curve, you should be able to obtain the first and last cells by doing the following:

first_cell_index <- which.min(fit$ord)
last_cell_index <- which.max(fit$ord)

However, I'm not sure this is so informative - if you're looking to get representative gene expression for the start and end of the curve, you should really be looking at the cluster of cells towards the start and end of the curve - for example, by comparing cells in clusters 1 and 3 that you showed in your first plot. A single cell will suffer from dropout and not tell you a lot about the cluster.


BERT Language Model and Gene Sequences - How Do I Relate Clusters of Sequences? I hope you'll indulge a question from a computer scientist with limited bioinformatics knowledge.   I've been working with the Google tool for language modeling called BERT.   It's generally regarded as state of the art when encoding language.   It allows you to say that the sentences "I went to the bank to cash my check." and "My wife is going to the Credit Union to pay the mortgage" are more similar than if compared to "he was fishing and fell off the bank."   In short, BERT takes sentences and encodes them so that related concepts are near each other and unrelated concepts are farther apart.

BERT is an amazing tool and I wondered if it could successfully model the language of the genome if it were encoded into "sentences" (I chose 7-mer words and sentences of length 64, arbitrarily).  Below is the result of encoding the full genomes of human, bonobo, gorilla, chimp, and orangutan.  There are about 270K "sentences" here.   You can see there are four major blobs, with indications of small tight clusters.   BERT has decided that the points indicated in red are somehow semantically related.  I think is is significant, or at the very least interesting.

My first thought was that clustering was trivial degenerate cases.  Maybe all the "AAAAAAA..." were put near each other, and all the "TTTTTTTTT...." etc.  I printed out samples from the clusters and visually that is not the case.  Also I've made histograms of the average Levenshtein distance in groups.  Within a tight group, the LD is statistically much lower (but not trivially low) as compared to unclustered regions (the big blob is all unclustered).

If there is a biological basis for the clustering, then new sequences could be similarly transformed, and classified according to cluster membership (and related biological function).  This exhausts my biological knowledge.  I've taken sequences from clusters and used BLAST to try and get some insight but I'm way out of my comfort zone with that.

My primary question is, given a set of sequences (in my case len=448), how do you figure out what they all have in common (given my hypothesis is that BERT says they are related based on proximity)?    And my backup question is, do you think this is interesting or a waste of time?

The UMAP results look very pretty :)

If I were to try to make sense of this blob, I'd try to do the following.

Pull out Gene Ontology (GO) terms for a list of candidate genes (you can try to subset Uniprot / Swissprot for some examples). Then I feed through some genes through your fragmented BERT kmer approach to see where they land. I'd be particularly interested in seeing where transmembrane proteins and enzymes land here.
I'd try to color this by introns vs exons (you should be able to extract this sort of information from your genome files). There is a good chance that the difference between introns / exons could explain some of these differences.
It maybe worth looking into genes that only encode RNAs vs genes that encode proteins.


Interpreting MultiQC plot of GATK BaseRecalibration data? I'm using GATK4 to build a somatic variant calling pipeline. The pipeline uses MultiQC to aggregate quality control data, and one of the QC measures reported is base quality score recalibration from GATK4's BaseRecalibraton. The description in MultiQC says

  This plot shows the distribution of base quality scores in each sample before and after base quality score recalibration (BQSR). Applying BQSR should broaden the distribution of base quality scores.

The black line is for my tumor sample and blue line is my normal sample. I don't understand how the before and after distributions are both shown in the plot as I only see 1 distribution per sample. Am I missing something, or should there be a "Post-calibration Count" plot that isn't shown? If nothing is missing, how should I be interpreting the plot?

Turns out two recalibration tables are needed to show pre-recalibration and post-recalibration plots. I previously only had generated the pre-calibration table by running BaseRecalibrator on the BAM that hadn't undergone ApplyBQSR to adjust base scores. Here are the steps to show pre and post recalibration plots in MultiQC:

Run BaseRecalibrator on your initial input BAM. This will create a recalibration table with the original base scores.
Run ApplyBQSR with the recalibration table from step 1 to create a new BAM with the adjusted base quality scores.
Run BaseRecalibrator on the adjusted BAM generated in step 2 to create a new recalibration table with the adjusted base quality scores.
Use both recalibration tables as inputs to MultiQC. This should create both plots.


Network comparison of single cells (from sequencing data) What methods can we use to compare networks (PPI, gene regulatory etc) created from single-cell gene expression (or proteomics) data?

There are methods that construct networks by comparing two conditions (e.g. from (gene) co-expression, or differential gene correlation (DGCA, Differential Gene Correlation Analysis)), and infer properties based on these networks.

What I'm interested here is the comparison of networks constructed for each single cell. The networks can be constructed solely from the data, or  using known interactions (edges) with the help of a database. These comparisons could be used for understanding or classifying cells.

This will probably require a network construction method (which may be hard because we do not have replicates for single cells), some form of network similarity measure / network alignment method. This is a very hard task in general, but here we have 1-to-1 correspondence for each node between the cells (networks), which should help.

Similar discussions at

Biostars: Question: How Do You Align Or Compare Two Networks
Content-Wise?
Researchgate: How can I measure similarity between two networks?

One related paper for single-cells is Chan et al., 2017, which uses groups of cells (created by other methods, e.g. clustering) to construct networks, and manually compares them (common edges / Venn diagram).

You can't build a network of a single cell only with the expression of a single cell. You either need previous known interactions or pathways or you need to use several cells/samples.
If you use previous known information, you can use pathway information, otherwise you can group some cells and use something along the lines of WGCNA to find a scale-free network.

In the first case one can add the expression upon it and see what is happening, what edges are kept and which ones not...

Once you have the network for each cell you can compare the metrics for each cell. Some metrics to look at are the centrality of each gene and the betweneess.
In a more general way, you can look if it follows a scale-free topology or other type of network structures, or how many nodes are present, how many edges are present...