MiXCR: only create a single export file for all clonotypes I am using following command from MiXCR to both align, assemble and export my input files:

mixcr analyze amplicon \
    -s hsa \
    --starting-material rna \
    --5-end no-v-primers \
    --3-end c-primers \
    --adapters adapters-present \
    --export "-aaFeature CDR1 -nFeature CDR1 -aaFeature CDR2 \
              -nFeature CDR2 -aaFeature CDR3 -nFeature CDR3 \
              -aaFeature VGene -cGenes -count" \
    --assemble "-OassemblingFeatures=[CDR1,CDR2,CDR3]" \
    INPUT_FILE.fasta analysis_name

Unfortunately, this creates multiple output file (analysis_name.clonotypes.ALL.txt, analysis_name.clonotypes.IGH.txt, etc.), which stands in contrast to the output of the mixcr exportAlignments command.

How can I export a single txt file instead, while still using the mixcr analyze command?

This is from the MiXCR documentation:

Under the hood analyze amplicon is equivalent to the execution of the following
MiXCR actions options:
> mixcr align -s hs -p kAligner2 input_R1.fastq input_R2.fastq alignments.vdjca
... Building alignments
> mixcr assemble alignments.vdjca clones.clns
... Assembling clones
> mixcr exportClones --chains IGH clones.clns clones.txt
... Exporting clones to tab-delimited file

So when using analyze amplicon, you are using mixcr exportClones under the hood and not exportAlignments, hence the difference in output.

I believe what you are experiencing results from a bug, please see this issue. If this is the case, all these files should be the same according to the bug report. Fortunately the developers have already addressed the bug in version 3.0.3.


How to convert fastq to fast5 fast5 is a variant of HDF5 the native format in which raw data from Oxford Nanopore MinION are provided. You can easily extract the reads in fast5 format into a standard fastq format, using for example poretools.

Say I have aligned these reads in fastq format to an external reference genome, resulting in a SAM file. Say I have then taken a subset of the  SAM file, according to the bitwise flag, to include only the reads that map to the reference. With the read ID, I can then grep them out from the file containing the reads in fastq format, generating a subset file in fastq format containing only the IDs that have mapped to the reference.

Now my question is, can we subset reads from the fast5 archive according to the list of mapping reads as taken from the file with reads in fastq format? This is for educational purposes, so that we have a smaller starting archive, and the fast5 -> fastq extraction takes less cpu time.

NOTICE:
I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.

This is something I have been meaning to get around to for a while, so thanks for the prompt.

I have created a python program called fast5seek to do what (I think) you're after.

As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.

The docs on the GitHub repo have all the info, but for those reading along at home

pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

What it does is read in <in.fastq|in.bam|in.sam> and extract the read id from each header. It then goes through all the fast5 files under /path/to/fast5s and checks whether their read id is in the set of read ids from <in.fastq|in.bam|in.sam>. If it is, the path to the file is written to it's own line in out.txt.

If no output (-o) is given, it will write the output to stdout.

So if you wanted to pipe these paths into another program, you could do something like

mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to subset_dir/.


Unable to use phangorn::phyDat I am following an example from ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.
He is doing:

I am doing pretty much the same:
x=c('garras','garras','garras', 'tejido','tejido','tejido')
y <- phyDat(matrix(x), "USER", levels = unique(x))

but I get:

Error in phyDat.default(data, levels = levels, return.index =
return.index, : data object must contain taxa names

Screenshot:

What am I doing wrong? how come the example from the book works...

I cannot access the book, it might be an older version or there are some steps in between. The screen shot doesn't help.
If you look at the vignette and the error, what you need to provide is a matrix of characters with defined row names, for example:
set.seed(111)
x = sample(letters,100,replace=TRUE)
phyDat(matrix(x),"USER",levels=unique(x))

Error in phyDat.default(data, levels = levels, return.index = return.index,  :
  data object must contain taxa names

mat = matrix(x,10,10,dimnames=list(paste0("t",1:10),NULL))
phyDat(mat,"USER",levels=unique(x))
10 sequences with 10 character and 10 different site patterns.
The states are n t s o y e c h z q m j d o h n z y x p q g u o a i h y c g t d f y t d f e u e f d i j r z x w u m a j p j g d t q y v o x p m s o z x z a y j z r u s a y k l f g l f r i l r f m t t e f x v f g q j


Gene expression signature (AUC) from gene expression data with 48 +48 samples Does anyone have a good method to find a signature where you combine the expression from multiple genes to predict a specific condition.
Assume you have performed sequencing of 48 cancer samples and 48 normal sample. Then find a combination of genes that predicts the cancer samples, using for instance AUC.

I have a method for doing this that uses a MCMC algorithm originally designed for population genetics (via the structure program) to generate a statistic between 0 and 1 indicating genetic risk. This needs the use of reference populations for cases and controls that the query individuals are added to for classification. I use T1D as an example. My methods paper is not peer-reviewed, and probably doesn't go into sufficient detail for what you are looking for, but at least indicates one way this can be done:

https://www.biorxiv.org/content/early/2017/01/31/104497

Based on the work I've done with this, using 48 cases and 48 controls with a very specifically defined type of cancer should be sufficient to create a generalised signature from human populations. If a generic cancer profile is needed, the populations are not representative of a general human population, or there are some cases in the control group, then things get trickier and require more individuals. On the other hand, 1500 per group is probably overkill.


How to add OVERLAP_BP field to VEP output? I'm using the Ensembl Variant Effect Predictor (VEP) to annotate VCF files, and I am trying to add some additional fields to the output, among these is OVERLAP_BP, which is listed here in available fields:

  OVERLAP_BP - Number of base pairs overlapping with the corresponding variation feature

It's clear that some (many? all?) of these additional fields are only enabled by adding extra options to the VEP command (as listed on the RHS of the options table).

For example, the ZYG field is made available by adding the --individual flag to the VEP call.

But OVERLAP_BP doesn’t appear in the output fields column on that page, and when I add it to the output using --fields, it is consistently empty.

What options are required to make the OVERLAP_BP output field available to me in VEP?

I also opened this as a Github issue, and received the following response:

  It seems there are some inconsistencies on our end. I'm looking into it.
  Until the fix, using --fields OverlapBP,OverlapPC works with any output format, while --fields OVERLAP_BP,OVERLAP_PC works only with default output format.

  A small note that this overlap data is available for structural variants only.


What is the significance of enrichment factor regarding machine learning methods? I'm fairly new to the field of bioinformatics and ran into a question while reading a paper I found on bioRxiv. The overall setting of the paper is using multi-task deep neural networks for kinase classification. The paper that I'm referring to is Kinome-Wide Activity Classification of Small Molecules by Deep Learning (Allen et al., 2019) for anyone who's curious.

I'll get straight into the specific figure that I'm having trouble understanding.

I know that questions like these are a bit discouraged on SE, but I'm having trouble on how to start regarding the interpretation of this figure.

Based on my understanding, the point seems to be that multitaks DNN's have the highest overall EF/EFmax ratio amounts, but what does this ratio signify?

I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is a concept "to measure the added value of a search tool over another one or over the homogeneous distribution in the genome population." Does this mean that multitask DNN's are "better" than the other models, due to its higher number of EF/EFmax points?

But if we look at the curves on the right, it seems to say that for fewer number of compounds, the ratio for multitask DNN's is lower...

I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is...

Where possible, it's best to use explanations provided in the original source. This is because definitions could be different, depending on who is writing. Here's some descriptive information from the paper that may be more informative:

  The enrichment factor (EF) is another metric we use and is very
  popular for use in evaluation of virtual screening performance. EF
  denotes the quotient of true actives among a subset of predicted
  actives and the overall fraction of actives.... To evaluate the
  maximum achievable enrichment factor (EFmax), the maximum number of
  actives among the percentage of tested compounds is divided by the
  fraction of actives in the entire dataset. To more consistently
  quantify enrichment at 0.1% of all compounds tested, we report the
  ratio of EF/ EFmax. This normalized metric is useful because EF values
  are not directly comparable across different datasets due to the
  maximum possible EF being constrained by the ratio of total to active
  compounds and the evaluated fraction as shown in the equation above.

So for the purpose of this paper, they're saying that this ratio is an appropriately normalised substitute for a plain enrichment factor. High enrichment factor == good; low enrichment factor == bad.

Given the information you've provided, I'd agree that the curves in this figure appear to suggest that the multitask DNN performs worse with smaller numbers of active compounds (where Naive Bayes or Logistic regression would be a better choice). It's hard to tell from the graph (and legend) what the curves represent, and the clipping at EF/EFmax = 0 doesn't help.

The curves have a strange shape, suggesting the number of active compounds is dependent on the EF/EFmax ratio, rather than the other way round (which would make more sense to me). If I consider the curves to be something like a smoothed mean of the active compound number, then any place where there are fewer dots (i.e. low EF/EFmax for Multitask DNN) will have more associated error. For Multitask DNN, I see a large gap between 0.0 and 0.2, with lots of outliers on the 0.0 line (for all methods), so presumably the 0.0 outliers are skewing the curves.

Given that this is a preprint, the generation of the curves isn't specified, and the other figures (especially Figure 3) seem to tell a different story (i.e. that multi-task DNN performs better at all scales), I'd recommend contacting the authors to ask / clarify what the curves represent, or ignoring the curves entirely.


How to retrieve the yeast locus tag from uniprot via sparql? Let's say I want to retrieve the yeast-specific identifier for a certain protein, in the example below it will be P00330 which I would like to link to YOL086C.

When I go to uniprot's sparql UI and type

PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT *
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?db up:database <http://purl.uniprot.org/database/SGD> .

}

it returns

protein db
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446

How would I now go on and find the YOL086C? When I add the line

?db ?p ?o

to the query above, I don't see anything that would directly lead me to the identifier.

I also tried to access another database

?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

but this returns all kind of databases while I am interested only in the FungiDB. How could I select this one?

How could I link the two IDs to each other?

Ok, figured it out by looking at the scheme uniprot provides:

PREFIX up:<http://purl.uniprot.org/core/>

SELECT ?protein ?yeast_id
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein up:encodedBy ?gene .
    ?gene up:locusName ?yeast_id
}

So one starts at the accession number in the middle and then uses the path to the actual locusName.

Also the path via the database works fine, I just did not use the seeAlso correctly (protein was only bound to db but not to db2, that's why it returned all entries):

PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT ?protein ?db ?db2
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?protein rdfs:seeAlso ?db2 .
    ?db up:database <http://purl.uniprot.org/database/SGD> .
    ?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

}

This will return the desired output (one would still have to process the db2 entry, of course):

protein db  db2
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446  http://purl.uniprot.org/eupathdb/FungiDB:YOL086C


How to extract all peptide sequences for a UniProt protein entry with an API? I know how to generally pull down information for a UniProt entry using the REST API, for example:
https://www.uniprot.org/uniprot/?query=id:Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab
I'm looking at this entry as an example. I'm able to pull down most of the attributes I want with a variation of the API query above, however, I'm stuck on one part.
For all human proteins, I want to pull down the peptide sequences that come from that protein.
In the example above, it would be the sequences of the peptides associated with this section:

Can this be done with the API or do I have to come with a scraping method for all human proteins where I (1) get the protein sequence of each entry, (2) get this table, (3) calculate where the peptides start and end from the sequence in step 1?
I feel like this alternative will be much slower, so I want to check I'm not missing anything obvious.

If downloading the complete database is not an option, you may able to use the instructions from this help page (based on the gff format for any given query):
https://www.uniprot.org/help/download_subsequences


Metagenomics: Identifying most common sequences I am working on a project and used the following command:
vsearch --derep_fulllength filtered_merged.fa -sizeout -relabel Uniq -output dereplicated_filtered_merged.fa

and got the following output:
87373926 nt in 203453 seqs, min 310, max 480, avg 352
Sorting 100%
10981 unique sequences, avg cluster 2.0, median 1, max 1287
Writing output file 100%

The output had provided me with the data that 10981 unique sequences have been identified. But I cant seem to identify how many reads of the most common sequence were present in the input data.
Any suggestions will be kindly appreciated!

According to the VSEARCH docs, since you have specified --sizeout your abundances have been written into the FASTA headers:

--sizeout
Take into account the abundance annotations present in the input fasta file (search for the pattern ’[>;]size=integer[;]’ in sequence headers). That option is active by default when rereplicating.
Add abundance annotations to the output fasta file (add the pattern ’;size=integer;’ to sequence headers). If --sizein is specified, each unique sequence receives a new abun- dance value corresponding to its total abundance (sum of the abundances of its occur- rences). If --sizein is not specified, input abundances are set to 1, and each unique sequence receives a new abundance value corresponding to its number of occurrences in the input file.


extracting last few lines of a growing bgzip compressed .vcf.gz file I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.

The merge process has been running for a couple of days, and I want to estimate how much longer it will take. The file is a growing .vcf.gz that is being bgzip -c compressed from a Linux pipe:

bcftools merge [...] | some processing and filtering | bgzip -c > final.file.vcf.gz

How can I know at which point is the file currently, to predict how much longer it will take to complete the merge?

Is there a way to do something like:

tail --with-some-magic final.file.vcf.gz | bgzip -d -c

So that I can read at which chromosome is currently on?

Is there a way to guesstimate the final size of the merge file and compare from the sum of sizes of the original files? E.g. the sum of compressed sizes originally (before filtering and munging) is 20G:

find $PWD -name "ALL.chr*.vcf.gz" | xargs ls -l | awk '{print $5}' | csvtk summary -H -f 1:sum | numfmt --to-si

This gives 20G, but the merge file is at 27G and still running.

Any ideas?

27GB is not huge. You can zcat the file and see the last few lines. It will take ~10 minutes. It is theoretically possible to peek a bgzf file, but no such tools are available.

By the way,

  I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.

Probably you should use bcftools concat, not bcftools merge. merge combines VCFs from different samples. concat concatenates VCFs for individual chromosomes; each input VCF contains the same set of samples.


how to repeat records in fastq n times efficiently? How to iterate/repeat a record n times in a fastq file using bioawk? I wrote a python code using biopython, but it is very very slow. So, I am wondering if I can get some help by using bioawk. Thank you so much!
for example, in the file.fq file, I'd like to repeat read_id1 for 100times, read_id2 for 100times, ..., and write them to a new file.

Not sure why you would want to do this, but if the aim is to repeat all records in the fastq/fasta files n times, the most efficient way I can think of is just to concatenate the orginal file n times, similar to this thread.
i.e adapting the accepted answer would look something like this (I have not tested this):
for i in {1..n};do cat file.fq >> file_ntimes.fq ; done

Update based on comments:
To replicate each fastq entry back to back, maybe this hacky unix command can be adjusted (say for n=4):
n=10
paste - - - - <file.fq | awk '{for(i=0;i<$n;i++)print}' | tr '\t' '\n' > file_10times.fq

The idea is to paste blocks of 4 lines (i.e. 1 fastq entry) into one line and replicate this line n times. tr is reversing the paste action (since the separation by paste is a tab character). This won't work if fastq file for some reason has tabs in it which I don't expect but didn't test either.
Note, for fasta file this would be paste - - rather than paste - - - -
I haven't tested this, but if it works for you I could polish it up if needed.
Update2
Based on Xioafeis comment this is solution seems to work as well:
bioawk -cfastx '{for(i=1; i<=10; i++) {print "@"$name"\n"$seq"\n+\n"$qual}}' file.fq > file_10times.fq


Intersection of two genomic ranges to keep metadata I am trying to find intersection of two genomic ranges (gr1 and gr2) and keep metadata from one of them

gr1
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

gr2
chrI [50, 150] + |

intersect(gr1, gr2) will remove metadata, where
subsetByOverlaps(gr1, gr2) return:

chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

where metadata is saved, but intersection is bigger than real (the same as in initial gr1 start/end)

How can I get an intersection range and keep metadata?
to get the range:

chrI [50, 100] * | 0.1
chrI [101, 150] * | 0.2

Presuming you don't have cases of multiple ovelaps

> library(GenomicRanges)
> gr1 = GRanges(c("chrI", "chrI"), IRanges(start=c(1, 101), width=c(100, 100)), mcols=data.frame(foo=c(0.1, 0.2)))
> gr2 = GRanges(c("chrI"), IRanges(start=c(50), width=c(101)))
> o = findOverlaps(gr1, gr2)
> grl1 = split(gr1[queryHits(o)], 1:length(o)) # You can't mendoapply on a GRanges object
> grl2 = split(gr2[subjectHits(o)], 1:length(o))
> foo = function(x, y) {
+     rv = x
+     start(rv) = max(start(x), start(y))
+     end(rv) = min(end(x), end(y))
+     return(rv)
+ }
> unlist(mendoapply(foo, grl1, y=grl2))

The output is:

GRanges object with 2 ranges and 1 metadata column:
    seqnames     ranges strand | mcols.foo
       <Rle>  <IRanges>  <Rle> | <numeric>
  1     chrI [ 50, 100]      * |       0.1
  2     chrI [101, 150]      * |       0.2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

It'd be more memory efficient to iterate over o, but this gives you the idea.

In cases of multiple overlaps, one could use intersect(), then findOverlaps() with that and then continue on with the procedure demonstrated above (or iterate over the overlaps).


How to search a specific sequence in BAM files for 10X experiment I have to search a specific sequence in a set of cells (> 1k cells) from a single-cell experiment done with 10X genomics.

As input file I have a single bam file, and 24 fastqs, therefore each file contains information about several cells. How can I search for a specific sequence and underpin which single cells have it?

I should be able to search for a specific sequence in a bam using samtools, etc. But the issue here is how to demultiplex the data from different cells to understand which cell provide the sequence.

The bam file has bam tags which say what reads belong to what cells.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam


Is there any difference between SNPs 'AG' and 'GA' in association analyses? I am using plink to check the association between a phenotype and the SNPs of a gene. plink says the phenotype is significantly associated with a SNP on that gene, and when I check the SNP alleles, for 6 out of 49 samples, the nucleotides are 'GA', and for the rest, they are 'AG', so only the order is different. I am not a biologist, but my basic genetics knowledge says that those two are semantically the same. However, obviously, plink does not agree with me and it associates this SNP with a small p-value since 'GA' samples have a higher value for the phenotype than 'AG' samples. I am very new to SNP analysis, any help is appreciated.

Could you please show us the context in which this appears, as you seem to be interpreting this differently to Devon.

If it's appearing as you say, GA and AG, then Devon is right, this usually means that the sequence goes ###AG### or ###GA###, which are two very different sequences.

If, however, as you're implying, the sequences are actually ###G### or ###A###, then semantically, the two ways of writing it have the same outcome. Normally, however, this is written as A/G or G/A.

Alternatively, you're looking at genotypes, which would usually appear as A|G and G|A. In this case, both individuals are heterozygotes, and in terms of the individual SNP they are the same. However the side of the pipe is relevant as this tells you which homologous chromosome it is on. This allows you to calculate LD with respect to neighbouring SNPs, and to determine the actual sequences of the proteins in that individual.


What would be a reasonable healthy data set/control group to compare with a cell line from TCGA? Assume I want to run a simple machine learning of healthy/non healthy classification, how can I use TCGA where it has only data of people with cancer. What do usually people do for example for study of gene expressions?

In the TCGA datasets there is a variable (type) containing the information if the sample was from the tumor or from an adjacent region, which is usually considered as healthy.

As this samples come from the same patient and the same tissues they are used in the analysis for the comparisons and classifications.


Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output? After alignment using BWA, I have removed the dupliment using the samtools(Version: 1.9).

My procedure is as follows:

bwa mem -k 32 -M ref.fa read1 read2 > out.sam
samtools view -@ 0 -b -T ref.fa -o out.bam in.sam
samtools sort -n -o out.nameSrt.bam in.bam
samtools fixmate -r -m in.nameSrt.bam out.fixmate.bam
samtools sort -o out.fixmate.sort.bam in.fixmate.bam
samtools markdup -r -S -s in.fixmate.sort.bam out.markdup.bam
samtools flagstat in.markdup.bam > out.markdup.flagstat

The flagstat output result is as follows:

21611397 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
21611397 + 0 mapped (100.00% : N/A)
21422330 + 0 paired in sequencing
10711165 + 0 read1
10711165 + 0 read2
19797684 + 0 properly paired (92.42% : N/A)
21422330 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
1306000 + 0 with mate mapped to a different chr
727043 + 0 with mate mapped to a different chr (mapQ>=5)

Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
Is there anything wrong with my procedure?

Bwa-mem may produce chimeric alignment, where different parts of a read are mapped to distinct loci. flagstat counts them as two reads.


Aligning sequence and comparing it against primer I am looking to show how a primer is consistent among some genomic data. I have a primer of about 23bp and looking to compare it to about 5000 genomic sequences of 10kb. I am unsure how to format it the way I need to, since I am just trying to show the nucleotide differences of the primer to the genetic sequences. Here is what I need:

>     1. Cut out the area that my primer is located and about 20bp down each end.
>     2. Show only the bases that are different from my primer in my analysis.
>      ex: Primer: -----------ATGTGGAAGCAAATATCAAATGA---------
>          Gene:   ATGACCATACG----C--------------T---ATCGTAGGG
>                  ATGAGCATACC-----A----T--------T---TTCGAACGC

The data I am using is all dengue sequences (all serotypes) and the primer with the following code: ATGTGGAAGCAAATATCAAATGA.

I am using Biostrings, msa, and seqinr. I use ncbi to get the genetic sequences and using FASTA files.

Thanks!

Edit: Solution:
I used MAFFT to align the data, working very efficiently and put it into a fasta file. Then I cut the nucleotides I did not need by:
1. Finding the part i wanted and cutting 10 nucleotides before it and after it.
2. Putting only the [start:end] into a new list, which I will then manipulate.

However, I am still having difficulty finding a way of showing my data nicely because I can't find a suitable package. Any advice?

Edit:
Hi all,
Just wanted to show the code I made to solve the problem.

    #So the we will be using a dataframe that has the sequences as a string
    #Then I created a temporary dataframe to store each sequence, with the primer at the top

    tempDF <- data.frame(matrix(primer1, nrow = 1, ncol = 1),
    stringsAsFactors = FALSE)

    for(i in sequenceAsString){
      #Temporary vector for each sequence after compared to primer.

      temporary1 <- ""
      temporary1 <- paste(temporary1, substring(i,1,11), sep = "")

      for(j in 12:(nchar(i)-20)){
        if(substring(i,j,j) == substring(toString(reversePrimer1), j-12,j-12)){
          temporary1 <- paste(temporary1, "-", sep = "")
        }else{
          temporary1 <- paste(temporary1, substring(i,j,j), sep = "")
        }
      }
      temporary1 <- paste(temporary1, substring(i,nchar(i)-19,nchar(i)), sep = "")
      tempDF <- rbind(tempDF,temporary1)
    }

So, I essentially take the sequence I chose apart and compare each character to that of the reverseprimer. I then replace the ones that are the same with an '-'

Thanks everyone for their help!

If I understand you correctly, after your most recent edit, you have 5,000 sequences that are 63 bp long and in each you simply want to replace the nucleotides from 21 through 43 with - unless that nucleotide does not match the primer. Since - is normally used to represent gaps, I recommend using . instead for matching nucleotides.

I know you're looking for an R solution, but I don't usually use R for string manipulations and don't have time at the moment to figure it out. So, assuming this might be urgent for you, below is a Python solution. This uses your example data which shows 11 nucleotides before the primer and 9 after it (so you'll need to adjust it to work with your actual case in which 20 nucleotides are before and after the primer). In Python, strings are arrays of characters, by the way.

Here are your example sequences are in the file seqs.txt:

ATGACCATACGATGTCGAAGCAAATATCAATTGAATCGTAGGG
ATGAGCATACCATGTGAAAGCTAATATCAATTGATTCGAACGC

Here is the code that gives you the expected output:

#!/bin/env python
import sys
primer = 'ATGTGGAAGCAAATATCAAATGA'

print('-' * 11 + primer + '-' * 9 )
with open('seqs.txt','r') as f:
    for gene in f:
        sys.stdout.write(gene[0:11])
        for gene_index in range(11,34):
            primer_index = gene_index - 11
            if gene[gene_index] == primer[primer_index]:
                sys.stdout.write('.')
            else:
                sys.stdout.write(gene[gene_index])
        print(gene[34:])

Running this gives you:

-----------ATGTGGAAGCAAATATCAAATGA---------
ATGACCATACG....C..............T...ATCGTAGGG
ATGAGCATACC.....A....T........T...TTCGAACGC


How to select only RNA with Hetero atoms from pdb file with python? I'm trying to separate RNA from protein in a complex protein/RNA PDB file and I want all RNA info with the hetero atoms in between the bases BUT without H20 etc. In short I want RNA part of pdb file without discontinuous lines.

I managed to separate RNA from protein with Bio PDB Select but it consider hetero atoms as amino acid when I use is_aa(residue). So hetero atoms wont appear in my "only RNA" file.

from Bio.PDB import *
from Bio.PDB import PDBParser, PDBIO, Select
import os

class ProtSelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == True else 0

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == False and residue.id[0] != "W" else 0

pdb = PDBParser().get_structure("2bh2", "pdb2bh2.ent")
io = PDBIO()
io.set_structure(pdb)
io.save("seqprotest.pdb", ProtSelect())
io.save("seqRNAtest.pdb", RNASelect())

EDIT: I managed to do it:

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

It's not pretty but it do the job.
The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H_' and ends by 'A,T,C,G,U'.
Just remove 'H_ MG' and you're good.

Solution I found: (c/p of body)

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

It's not pretty but it do the job. The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H_' and ends by 'A,T,C,G,U'. Just remove 'H_ MG' and you're good.


Mapping Information from SAM/BAM file I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM)

How many PacBio reads are mapped to at least one illumina reads
A Table showing each pacbio read, the number of illumina reads that
mapped, and total coverage (i.e. percentage of mapped regions)

If you sort and index your BAM file:

samtools idxstats. The third column is the number of Illumina reads mapped to each PacBio read. You can then pipe that to awk: samtools idxstats foo.bam | awk '{if($1 != "*" && $3 > 0) tot += 1}END{print tot}. That will print the number of PacBio reads with at least 1 read mapped to them.
You'll want samtools depth for that

samtools depth -aa foo.bam | awk '{if($1 != lread) {if(lread != "") {print lread, 100*lcov/llen}; lread = $1; lcov=0; llen=0}; llen += 1; if($3>0) {lcov +=1}}END{print lread, 100*lcov/llen}'

That won't give you the number of illumina reads per pacbio read, but you can get that from samtools idxstats (column 3) and join the results.

BTW, your usage of "coverage" is non-standard. You want "percentage coverage" rather than "coverage", which actually means "fold-coverage" (i.e., number of reads times their average length divided by the total genome/chromosome/etc. length).

If you instead want "fold-coverage" then the easiest way is to first know the length of your reads. Presuming you have 100 base reads (for paired-end reads, this would be the same, presuming the ends don't overlap):

samtools idxstats foo.bam | awk '{print $1, $3, 100*$3/$2}'

That will give you a table of PacBio read name, number of aligned reads, fold-coverage.


Adding an attribute to GFF3 file I failed to add Note=Gene description to mRNA attribute with the below code:

#!/usr/bin/python3
import click
import gffutils
import gffutils.gffwriter as gffwriter

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    #db = gffutils.create_db(gff3, dbfn='data/test.db', force=True, keep_order=True, merge_strategy='merge',
    #                        sort_attribute_values=True)

    db = gffutils.FeatureDB('data/test.db', keep_order=True)
    for gene in db.features_of_type('gene', order_by='start'):
        print(gene)
        for i in db.children(gene, featuretype='mRNA', order_by='start'):
            print(i)
            print(i.id)
            print(i.attributes['Note']=["Gene description"])
            db.update([i])

I got the following error:

File "/projects/bioinf-scripts/mod_braker_gff3.py", line 20
    print(i.attributes['Note']=["Gene description"])
         ^
SyntaxError: keyword can't be an expression

What did I miss?

Thank you in advance

You are printing the attribute before setting it, that's why it's complaining. So instead of this:

    print(i.attributes['Note']=["Gene description"])

Set it first and then print:

i.attributes['Note']=["Gene description"]
print(i)


gffread: GFaSeqGet errors on coordinate overhang Disclaimer: I had this issue posted on Tuxedo Tools users group and shared it on Twitter, but could not get an answer to this from the developers, nor find documentation of this issue online. So, I'll share my solution here below and see if some of you have a better solution to add.

Using gffread utility to convert transcripts from GTF as output of Cufflinks to FASTA, I get some GFaSeqGet errors, as follows (here is the full log file):

$ gffread transcripts.gtf -g ../genome.fa -w transcripts.fasta
Error (GFaSeqGet): end coordinate (32519) cannot be larger than sequence length 32512
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (6327) cannot be larger than sequence length 6326
Error (GFaSeqGet): end coordinate (4651) cannot be larger than sequence length 4648
Error (GFaSeqGet): subsequence cannot be larger than 3551
Error getting subseq for CUFF.29386.2 (1..3568)!

That's a bit weird since I used the same genome file as the one used to generate the alignment file required by Cufflinks, and the coordinates overhang does not make a lot of sense to me.

The alignment file (SAM) was generated using Hisat2:

$ hisat2
    --rna-strandness RF
    --dta-cufflinks
    -x ./hisat_idx
    -q
    -1 ../1A_1.fq.gz,../1B_1.fq.gz,../2A_1.fq.gz,../2B_1.fq.gz,../3A_1.fq.gz,../3B_1.fq.gz
    -2 ../1A_2.fq.gz,../1B_2.fq.gz,../2A_2.fq.gz,../2B_2.fq.gz,../3A_2.fq.gz,../3B_2.fq.gz
    -S out.sam

Then sorted:

$ cat out.sam | sort -k 3,3 -k 4,4n -T . > out_sorted.sam

My question is: are only the few sequences reported in the log affected? Or is it a more general problem? The answer seem to be the second one.

Out of gffread, alongside the log file, I actually got a FASTA file with no additional error messages during its generation (537M out of 404M GTF file) - but there are around 13k transcripts which are present in the GTF but seem to be missing from the FASTA:

$ cat transcripts.gtf | awk '$3=="transcript"{print $0}' | wc -l
231139
$ cat transcripts.fasta | grep '^>' | wc -l
218075

Do you know what would be a cause of this problem?

Basically, the source of the problem is issues with hard/soft clipped bases in short-reads RNAseq reads alignments to reference genome in post-transcriptome assembly steps.

Truncating the coordinates in the gtf file to match the reference is a solution. For this, I will share a perl solution (non-official) that will do the job.

perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf > transcripts_fixed.gtf

Where:

genome.fa.fai is the fasta index of the genome file
transcripts_to_fix.gtf is the transcript set out of the assembler with issues
transcripts_fixed.gtf is the new transcript file with the repaired end coordinates

Content of cufftrim.pl:

#!/usr/bin/perl
use strict;
use Getopt::Std;

my $usage = q/Usage:
 cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
umask 0002;
getopts('o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
if ($outfile) {
  open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
  select(OUTF);
  }
# --
my $fai=shift(@ARGV) || die ("$usage\nSamtools fasta index required!\n");
my $gtf=shift(@ARGV) || die ("$usage\nCufflinks' transcripts file required!\n");
my %slen;

open(FAI, $fai) || die("Error opening $fai\n");
open(GTF, $gtf) || die("Error opening $gtf\n");
while (<FAI>) {
 chomp;
 next if m/^#/;
 my @t=split(/\t/);
 die("Parsing error: $t[1] length for sequence $t[0] ?\n")
    unless $t[1]>10;
 $slen{$t[0]}=int($t[1]);
}
close(FAI);

# now process the GTF file
my ($tcount, $ecount); #trim counts
while (<GTF>) {
 my @t=split(/\t/);
 my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
 if ($t[4]>$l) {
   $t[4]=$l;
   $t[2] eq 'exon' ? ++$ecount : ++$tcount;
 }
 print join("\t", @t);
}
close(GTF);
print STDERR "$tcount transcripts, $ecount exons adjusted.\n";
# --
if ($outfile) {
 select(STDOUT);
 close(OUTF);
 }

Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution.


Biopython BioSQL error : AttributeError: 'DBSeq' object has no attribute '_data' Hi I am playing with Biopython and BioSQL,
I manage to create an MySQL database using Managing local biological databases with the BioSQL module tutorial:
and load 2909 fasta sequences in it through:
db = server["test_1"]

count = db.load(SeqIO.parse("Newtextfile-03.txt", "fasta"))
        print("Loaded %i records" % count)
        server.commit()

from this file: Newtextfile-03.txt
problem is trying just to print the database index list with:
for i in db:
    print(i)

works , while trying print element by element:
for i in db:
    print(i, db[i])

or with:
for key, record in db.items():
            print ("Key %r maps to a sequence record  %s" % (key, record))

gives an error:
2557 Traceback (most recent call last):
  File "./try_001.py", line 142, in <module>
    f()
  File "./try_001.py", line 120, in f
    print(i, db[i])
  File "......../lib/python3.8/site-packages/Bio/SeqRecord.py", line 652, in __str__
    lines.append(repr(self.seq))
  File "........//lib/python3.8/site-packages/Bio/Seq.py", line 109, in __repr__
    return f"{self.__class__.__name__}({self._data!r})"

AttributeError: 'DBSeq' object has no attribute '_data'

the error occurs for three records namely : [2164, 2522, 2557]
inspecting the Newtextfile-03.txt file results in finding fasta sequences of less than
100 aminoacids.
Is there any way to resolve the error ? It sounds strange to me that I can load the records into the database but not retrieving them. Is anything I am missing about the record load ? or similar ?
I am working with the database using:
 server = BioSeqDatabase.open_database(
    driver='MySQLdb' ,
    user=password.userz,
    password=password.pazzword,
    host="XXXXXXXX",
    db="bio-data-first")

db = server.new_database("test_1", description="Just for testing")
        server.commit()

using  driver="MySQLdb", wasn't able to have driver='mysql.connector' to work,
biopython is version 1.78

It's a bug in Biopython 1.78.
What those three sequences have in common is that they're less than 61 characters long.  When you print the record it implicitly calls __str__ for the record, which calls __repr__ for the Seq.  DBSeq gets its method from Seq, which looks like:
    def __repr__(self):
        """Return (truncated) representation of the sequence for debugging."""
        if len(self) > 60:
            # Shows the last three letters as it is often useful to see if
            # there is a stop codon at the end of a sequence.
            # Note total length is 54+3+3=60
            return f"{self.__class__.__name__}('{str(self[:54])}...{str(self[-3:])}')"
        else:
            return f"{self.__class__.__name__}({self._data!r})"

...and then it only tries to access the nonexistent _data for the short sequences.
So there's a mismatch between DBSeq's lack of a _data attribute and its use of Seq.__repr__.  I don't see a DBSeq class at all in the latest Biopython, so maybe this is fixed already.  But there's no tag for 1.79 yet that I can see and I can't get the very latest code to import successfully so I'm not sure.  In the meantime if you avoid letting repr(record) happen you should be able to work around it.


Files for paired end RNA sequencing I am looking at the videos at a DIY Transcriptomics course and the speaker mentions that to run Kallisto for read alignment with paired end sequencing, one would enter the following:

where the sample files at the end refer to the two different end reads (forward and reverse read files). However, browsing the NCBI for paired end data, I noticed that only 1 Fastq file is provided per replicate for a paired end procedure (here's an example). So I'm a bit confused about this.

The problem is that the SRA website can be a little confusing. View the same project on ENA where the organization is more human-readable: https://www.ebi.ac.uk/ena/browser/view/PRJNA574273


Getting weird plot with ggpubr package I have a data frame for two groups and 9 variables

> head(aa)
  treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
> class(aa)
[1] "data.frame"
> str(aa)
'data.frame':   1026 obs. of  3 variables:
 $ treatment: Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
 $ variable : Factor w/ 9 levels "Sig1","Sig2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value    : num  0.117 0.189 0.219 0.128 0.171 ...
> tail(aa)
     treatment variable value
1021      post     Sig9     0
1022      post     Sig9     0
1023      post     Sig9     0
1024      post     Sig9     0
1025      post     Sig9     0
1026      post     Sig9     0
>

I am using this code

library(ggpubr)
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> my_comparisons <- c("pre","post")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 50)
Warning message:
Computation failed in `stat_signif()`:
missing value where TRUE/FALSE needed
>

But I am getting this weird plot

While something like this is my goal

I really don't know where I am doing wrong

I used ggplot where I got something like this but I need statistics that is why I have stuck on this

I then used this code

> DF <- aa %>% pivot_longer(.,cols = c(Sig1,Sig2,Sig3, Sig4,Sig5,Sig6,Sig7,Sig8,Sig9), names_to = "var", values_to = "val")
> View(DF)
> library(ggplot2)
> library(ggpubr)
> ggplot(DF, aes(x = treatment, y = val, fill = treatment))+
+     geom_boxplot()+
+     stat_compare_means(comparisons = c("pre","post")+
+     stat_compare_means()+
+     facet_grid(var~., scales = "free")
+

But I got another weird plot

But changing the scale I got this but gives one statistics

While I need one stat for each variable between pre and post groups

This solution gives to statistics

> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment",facet.by = "variable")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y =2)

When I edited my code

> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> p + stat_compare_means(aes(group = variable))
>

But again one stat

The problem is the scale used: For the plot you called "weird" (first from the top), the scale is 50 and for the "ggplot only" (third from the top) the scale is 1.

You should play with the stat_compare_means(label.y = 50) bit, you can try setting the label.y parameter to 1.5 or 2.

EDIT BECAUSE OF THE ADDITIONAL REMARKS ON THE QUESTION:

Once again (and for sure the last time), I had to manually create example data as OP ignores repeated reminders on how to make a reproducible example with R.

library(data.table)
library(ggpubr)

text <- "
treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
1       post     Sig1 0.1173863
2       post     Sig1 0.1888243
3       post     Sig1 0.2191765
4       post     Sig1 0.1277735
5       post     Sig1 0.1711712
6       post     Sig1 0.0000000
1021      pre     Sig9     0
1022      pre     Sig9     0
1023      pre     Sig9     0
1024      pre     Sig9     0
1025      pre     Sig9     0
1026      pre     Sig9     0
1021      post     Sig9     1
1022      post     Sig9     1
1023      post     Sig9     1
1024      post     Sig9     1
1025      post     Sig9     1
1026      post     Sig9     1
"
my_data <-fread(text)

> my_data
      V1 treatment variable     value
 1:    1       pre     Sig1 0.1173863
 2:    2       pre     Sig1 0.1888243
 3:    3       pre     Sig1 0.2191765
 4:    4       pre     Sig1 0.1277735
 5:    5       pre     Sig1 0.1711712
 6:    6       pre     Sig1 0.0000000
 7:    1      post     Sig1 0.1173863
 8:    2      post     Sig1 0.1888243
 9:    3      post     Sig1 0.2191765
10:    4      post     Sig1 0.1277735
11:    5      post     Sig1 0.1711712
12:    6      post     Sig1 0.0000000
13: 1021       pre     Sig9 0.0000000
14: 1022       pre     Sig9 0.0000000
15: 1023       pre     Sig9 0.0000000
16: 1024       pre     Sig9 0.0000000
17: 1025       pre     Sig9 0.0000000
18: 1026       pre     Sig9 0.0000000
19: 1021      post     Sig9 1.0000000
20: 1022      post     Sig9 1.0000000
21: 1023      post     Sig9 1.0000000
22: 1024      post     Sig9 1.0000000
23: 1025      post     Sig9 1.0000000
24: 1026      post     Sig9 1.0000000
      V1 treatment variable     value

ggboxplot(my_data,
          x = "treatment",
          y = "value",
          color = "treatment",
          facet.by = "variable") +
  stat_compare_means(label = "p.format")

And this code generates:


Phasing problem using "samtools phase" so I'm having a problem using samtools phase. I'm trying to phase my .bam archives that where generated by aligning my samples with my reference genome. Following the documentation I tried the command:

samtools phase -A  my_referenceGenome.fasta   my_sample_toBePhased.bam

The -A flag stands for dropping reads with ambiguous phase. The error I get is: [W::sam_read1] Parse error at line 1. I generated .bam files with and without the SAM header but the error persists.

I feel like I am missing something basic but that isn't clearly stated anywere. Does the my_referenceGenome.fasta archive needs to be previously phased? If so, how can I do it?

Thank you so much!

This is what I see in my samtools phase (v1.10)

Usage:   samtools phase [options] <in.bam>

samtools phase only takes a bam file as an input. But you are specifying a .fasta file as an input. If you want to specify the reference genome, you should be using the argument "--reference" as described here

      --reference FILE
           Reference sequence FASTA FILE [null]


Is there any R package for biodbnet? I want to convert Ensembl ids to Gene_symbols. I used biomaRt package but all Ensembl ids are not seen in the results.

Examples of problematic IDs are: ENSG00000255414, ENSG00000274358, ENSG00000235565, ENSG00000278309. I gave them in bioDBnet, which seems to handle them.

I would like to know whether there is any R package for bioDBnet or not?
If so, could you please tell how to use it?

I couldn't find a package neither in Github nor in Bioconductor. Neither could I find any paper about a package linking to bioDBnet. So I don't believe there is any (R) package at the moment which links to biodbnet in order to retrieve the IDS.
However it has a couple of APIs so you could build one or use them to retrieve the information you want.


how do I predict bacterial small non-coding RNA for a specific mRNA? I am working on Vibrio parahaemolyticus. I have a gene that it is possible might be regulated by small non-coding RNAs. How can I predict possible sRNAs that target the transcript of this gene. Most software offers the opposite of what I want, i.e., they predict targets for small RNAs!
I'm trying to predict bacterial small non-coding RNA in general (Support Vector Machine) but I got more than 8 thousand sequences within intergenic regions alone. Is there a code for aligning these sequences to my transcript?

Answer from @llrs, converted from comment:
You could try to download RNAseq data from this bacteria and see if some of the RNA align to your sequence. You could try to blast your gene and look for the sequences that have low identity similarity (or classify them by type and then look what is each hit).


Save output to a specific folder and/or with a specific prefix in Cancer Genomics Cloud I would like to run the CNVkit workflow on Cancer Genomics Cloud.

However, the two workflows that I'd like to use:

CNVkit Reference, for building a reference from normal samples
CNVkit CNV Calling, for calling CNVs by comparing tumors to the reference

don't have parameters that can be set to define in which folder to save the results or which prefix to add to the results. This means that if you're generating multiple references or calling CNV for multiple tumor samples you will end up with results having the same name for each run (tumor).

I tried to modify the CWL workflow, online and through Rabix composer, but I was unable to make it work. I would like to use the published workflow to not code everything from the ground up.

Do you know how to update/use the current CWL workflow in order to avoid having the same filename/folder for all the results?

If the output is File, you can add this step https://github.com/Barski-lab/workflows/blob/master/tools/rename.cwl to rename the output file.


How I can reproduce this heat map I have single cell RNA-seq on 8 time points. I have clustered each time point to 2-3 clusters by seurat. I also have a list of differentially expressed genes between all my 1562 cells in 8 time points. I know how to plot this list of genes in heat map on each cluster of cells or a big heat map of my 1562 cells, but I have also sub-clusters for each cluster (time point).

Could you please help me to figure out how can I have a heat map with sub-clusters for each time point on which I am showing a list of differentially expressed genes?

Something like  the picture below where each day has been clustered to 2-3 sub-population then they plotted common differentially expressed genes between day points

Edit

Lets say 500 genes commonly are being differentially expressed in 8 time point. time point 1 has 3 sub-clusters of cells and time points 2 has 2 sub-clusters of cells. time point 3 and 4 both have 2 sub-clusters of cells. I want to heat map these 500 genes in cells in time point 1 (3 sub-clusters), time point 2 (2 sub-clusters), time point 23 (2 sub-clusters) and time point 4 (2 sub-clusters) in one unified heat map on which cells in each subclusteres (a,b and c) are in column and genes in rows but I don't have any idea how to do that. column annotation function in ComplexHeatmap likely works in this regard But I don't know how to define cells inside each  sub cluster in different colours

I can heat map of 500 genes only on time points easily by seurat but without underlying sub clusters

> DoHeatmap(object = SubsetData(object = object, max.cells.per.ident =
> 100), genes.use = 50 genes, slim.col.label = TRUE, group.label.rot =
> TRUE)

Figure 2 in this article.

In case it is still helpful as the post is rather old, the code below would generate a heatmap with annotations thanks to the ComplexHeatmap package.

But before that I would like to stress out the importance of example data, I bet a lot of experts could not help you just because of they did not have single cell data at their disposal.

I am using some scRNA-seq data stored within a Seurat object. The expression data come from the @scale.data slot and the annotations are from the @meta.data slot. I am using a few general cell type markers as genes and 3 columns from the metadata slot (cell type, cell subtype, sample region). When using the ComplexHeatmap package, the trick is to have the same order in the annotation and actual data tables and this is how you should control what (in your case single cells) is plotted in what order.

library(Seurat)
library(dplyr)
library(ComplexHeatmap)
library(circlize) # for the colorRamp2() function

seurat_object <- readRDS(...)

# Annotations

meta_data <- seurat_object@meta.data

# order of annotations/colors are defined here
ordered_meta_data <- meta_data[order(meta_data$cell_type), ]

# OPTIONAL: YOU CAN PICK COLORS FOR EACH LEVEL OF ANNOTATION
# HERE I PROVIDE COLORS FOR ONLY TWO OF THE THREE LEVELS, COLORS FOR THE
# REMAINING LEVEL IS TAKEN CARE OF BY THE PACKAGE
annotation_colors <- list("cell_type"= c("Epithelial" = "red",
                                         "Immune" = "orange",
                                         "Other" = "gray",
                                         "Stroma" = "green"),
                          "sample_region" = c("B" = "blue",
                                              "C" = "yellow",
                                              "N" = "maroon",
                                              "T" = "darkblue"))
ha = HeatmapAnnotation(df = ordered_meta_data,
                       show_annotation_name = TRUE,
                       col = annotation_colors)

# Expression data

genes_to_use <- c("EPCAM",
                  "CD3D",
                  "NCR1",
                  "CD79A",
                  "LYZ",
                  "COL1A1",
                  "CLDN5",
                  "S100B",
                  "KIT")

seurat_object <- ScaleData(seurat_object,
                           genes.use = genes_to_use)

my_data <- seurat_object@scale.data

# COLUMN ORDER OF THE EXPRESSION DATA SHOULD MATCH THE ROW ORDER OF THE
# ANNOTATION TABLE
my_data <- my_data[, rownames(ordered_meta_data)]

# Heatmap

col_fun = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
png(filename = "example_heatmap.png",
    width = 1000,
    height = 1000)
Heatmap(
  my_data,
  col = col_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  column_order = NULL,
  show_row_dend = FALSE,
  show_column_dend = FALSE,
  show_row_names = TRUE,
  show_column_names = FALSE,
  use_raster = TRUE,
  raster_device = c("png"),
  bottom_annotation = NULL,
  top_annotation = ha
)
dev.off()


error in heatmap using R I am learning about generating heatmaps using R. For this example I am using randomly generated data to test the process, but I encountered what appears to be an error or perhaps there is something am missing. Specifically, the values of Npy2r in groups A & B are the same, yet the heat map colors them differently. I am trying to understand why. Am including the R code I am using to generate the heat map

heatmap.2(x, Rowv=F, Colv = F, scale="col", col=greenred(3), key=TRUE, symkey=FALSE,
      density.info="none",cexRow=1,cexCol=1,margins=c(6,11),  trace="none",srtCol=45)

You have scale="col" in the code. What you are plotting is the z-score calculated based on the value distribution per column. Try changing it to "row" or "none".


Secretome and membrane receptor profiles I am looking the secretome profile and the membrane receptor profile for a given cell type.

In my specific case, this should be the secretome and outer membrane receptor profiles of dorsal root ganglion.

What I've done in the past is taken proteomics datasets or RNAseq datasets and used this as a reference cell surface proteome or secretome.

I wasn't sure what the best way to filter secreted proteins or membrane receptors were, so I used the Human Protein Atlas and this turned out to be a big headache. A possible alternative method might be to use DAVID, but there must be a more efficient method via the command line.

My questions are these:

Is there a "consensus" database containing protein expression
profiles for different tissue types in human?
If not, what is the best way to filter out the secreted and transmembrane proteins by gene symbol?

I don't think there is a database with expression data for your particular cell type. I am afraid your "headache" approach doesn't sound so bad...

So get expression data e.g., from GEO for your cell type, and define which genes/proteins you consider as "expressed" (e.g., minimal RPKM value).

Then, you can use GO annotation to find your genes involved with membrane or secretion (or what ever). You can use biomaRt in R for this.

For example, you want to know all "membrane" genes.

> library(biomaRt)
>
> ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
>
> membrane.genes <- getBM(attributes=c('hgnc_symbol',
> 'ensembl_transcript_id', 'go_id'),
> filters = 'go', values = 'GO:0016020', mart = ensembl)

And now you want to know it only from genes that are expressed in your ganglion cells (with dummy example of two genes).

> gene.list <- c("KIR2DL3", "GPR1")
>
> membrane.genes.list <- membrane.genes[membrane.genes$hgnc_symbol ==
> gene.list,]
>
> membrane.only <- membrane.genes.list[membrane.genes.list$go_id ==
> 'GO:0016020',]


bigWig to bed for regions above/below threshold I'd like to threshold a bigWig file on a value and get all regions above or below that value in bed format. Ideally without too many file conversions. Are there any tools to do that?

You can do that with a bit of python:

#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.


No MQ tags in VCF files To call minority variants in my Mtb sequences I'm using a pipeline of

fastq ->
    BWA MEM for mapping ->
    picard MarkDuplicates ->
    GATK3 for indel realignment and base recalibration ->
    lofreq for variant calling (https://csb5.github.io/lofreq/)

However the resulting VCF files don't have an MQ tag for me to filter on. Where in this process would mapping quality information be generated or lost so that it doesn't end up in the VCF file?

As @ATpoint explains in his comment, the mapping quality (MQ) applies to read mappings and, hence, is contained in BAM files, but not VCFs. The MQ values carry important information that ideally the variant caller should use when determining the variant quality (QUAL in the VCF) of individual variant loci.
In the LoFreq documentation it says that

Choice of mapper: Bowtie and BWA-SW
BWA-SW assigns very low mapping qualities to mapped reads. LoFreq incorporates mapping quality by default into its variant calling model, which doesn’t help in this scenario and will reduce its sensitivity. Consider disabling the use of mapping quality for lofreq call with -N.

So LoFreq is already aware of the MQ values and takes them into consideration during variant calling.
Since you are using BWA-MEM (and not BWA-SW) you should probably not be concerned about the issue they mention and trust their quality filtering.
If you are sceptic about the quality of mappings, though, you could still filter them out directly after mapping, e.g. using samtools view -q 20 ...


Which identifier do I have to use when I want to add a ##Fasta section to a gff? I have a gff file like this

scaffold_1  source  mRNA    2987526 2992430 .   -   .   ID=protein_id_68892;Name="foobar";locus_tag=my_organism_68892;translation=MHTGDALEGSTGNVSILV...
scaffold_1  source  exon    2987526 2987805 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2987526 2987805 .   -   1   name "foobar"; proteinId 68892; exonNumber 5;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2988610 2988745 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2988610 2988745 .   -   2   name "foobar"; proteinId 68892; exonNumber 4;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2989157 2989645 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2989157 2989645 .   -   2   name "foobar"; proteinId 68892; exonNumber 3;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2990903 2991541 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2990903 2991541 .   -   2   name "foobar"; proteinId 68892; exonNumber 2;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2991587 2992430 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2991587 2992430 .   -   0   name "foobar"; proteinId 68892; exonNumber 1;locus_tag=my_organism_68892;Parent=protein_id_68892

The mRNAentry was not present at first and I added it in. I also have the protein fast file available and would now like to add the amino acid sequences to the gff. I originally added them in column 9 to the mRNAentry which is apparently incorrect. Here it is written that the sequences should go at the end of the file in the ##Fastasection. However, I struggle to figure out which identifier I should use to link it to my mRNA entries.

Do I have to use

##Fasta
>protein_id_68892
MHTGDALEGSTGNVSILV...

or does it work with the locus_tag

##Fasta
>my_organism_68892
MHTGDALEGSTGNVSILV...

or is it something completely else? Where and how do can I add the translation information?

In the link above (and also in the comments it is suggested) the identifier in the first column is used, in my case that would be scaffold_1. However, I do have a lot of entries with scaffold_1in my file so a link between scaffold_1and an amino acid sequence would not be unique as I would have

##Fasta
>scaffold_1
MHTGDALEGSTGNVSILV...
>scaffold_1
GDAHTEGSTGNVSILV...

. In the end, I want to create a genbank file (which does not exist yet) where I can link a CDS to an amino acid sequence.

Unless I'm mistaken, the only variant of GFF that explicitly permits Fasta entries is GFF3, and the data you've shown above is definitely not valid GFF3. The mRNA entry seems to correctly use the = symbol to separate keys from values in each attribute, but the remaining columns all use spaces and quotation marks (GTF style) instead. This is likely to cause all kinds of problems with software packages that take GFF data as input.

But to more directly address your question, there is no "right" way to connect feature entries to sequence entries in a GFF3 file. There are some conventions, but these are all very specific to a small set of tools and community of practice, and it is rare that these conventions are objectively enforced. Some options include:

Placing genomic nucleotide sequences in the ##FASTA section, and then using feature coordinates to pull out mature mRNA sequences and translate them in silico to derive the corresponding peptide sequence. In this case, the Fasta entries would be genomic sequence identifiers that correspond to the first field/column in each feature entry.
Placing mature mRNA sequences in the ##FASTA section, and then directly translating this sequence to derive the corresponding peptide sequence. In this case, you'll want the sequence IDs in the ##FASTA section to match the ID attribute of the corresponding mRNA feature.
Placing the peptide sequences directly in the ##FASTA section. In this case, you would either want to have protein features declared in the GFF3 file (uncommon) or have some sort of protein_id attribute associated with each mRNA feature, in which case you'd want the sequence IDs in the ##FASTA section to match this protein_id attribute.

There are probably other ways of doing it, too. Which way is "right"? The GFF3 specification has nothing to say about it. It depends on the tools you are using to do the subsequent processing and analysis. If you're lucky, those tools might even document what they expect in terms of input.

Let me just say, as a post-script, that I think placing sequences and feature annotations in the same file is a horrible idea. I've never liked this about GFF3.


VerifyBamID freemix score We all have used freemix score at some stage to check contamination or swaps in our sequencing experiment, but can anyone once for all explain how the score is calculated?

Here's my current understanding of the freemix metric output by verifyBamID. The freemix score is a measure of contamination that is calculated without knowing the genotypes of the individuals. In the paper it is equation (2), and the authors refer to it as the "sequence-only" method. However, in addition to the sequencing reads, it requires some knowledge of SNPs and their allele frequency from the population. They demonstrate in the paper that using the SNPs and allele frequency estimates from a different population actually reduces the estimated level of contamination.

As for intuition for what equation (2) is measuring, imagine what would happen to the number of heterozygous sites detected if DNA from two different individuals was sequenced together. The total number of heterozygous sites would increase because it would count all sites where the two individuals are homozygous for opposite alleles or at least one of the two individuals is heterozygous. This would lead to an increased fraction of heterozygous sites over that expected from Hardy-Weinberg Equilibrium (i.e. 2pq).

Here's the explanation of freemix directly from the author, Hyun Min Kang:

  The key idea of FREEMIX estimate is to use excessive heterozygosity to estimate the level of contamination. Especially for common SNPs, you will observe higher fraction of heterozygous alleles than 2*p*(1-p), and it turns out that you can quantify the contamination very well if you know the population allele frequency already. If you do not have accurate population allele frequency information, than it would be harder to estimate FREEMIX parameters using verifyBamID.

Additional sources:

Verifying Sample Identities - Implementation
Interpreting output files
verifyBamID Google Group


How to normalise scRNASeq data for differential expression analysis I wish to perform differential expression analysis for cluster-specific gene expression in single-cell data (with a tool such as MAST or SCDE).

I have data for 3 biological replicates. I performed analysis on 10x data with cellranger and have filtered out cells with low UMIs. I have corrected for batch effects and identified clusters by performing graph-based clustering on the UMAP projections.

As I understand it, the components from batch effect correction algorithms are only intended to be used for the purpose of plotting and clustering. They are not to be used for analysis of individual gene expression or differential expression. This has been stated unequivocally by the original authors of MNN and CCA batch correction algorithms.

I have derived a UMI count matrix for each sample which corrects for amplification bias. Differences between the cells and samples still persist in the raw data. These are obvious from the mean UMI counts to the UMAP of uncorrected data. The number of UMIs differs considerably between samples (batches) and individual cells within a sample (the total UMI per cell). The amount of data from each of the cells varies considerably. The read-depth between samples is also different. Therefore, I need a simple approach to correct for differences in read-depth, total UMI or sensitivity between cells.

Is it possible to normalise UMI counts between cells? Should UMI count data be treated differently to reads per transcript in bulk RNA-Seq? Would it be more appropriate to do so at the level of individual cells or for batches across samples?

I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do:

markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch')

(change "object" and "batch" accordingly)

See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's group for some more info on the LR test.


Hemoglobin subunits genes in scRNA-seq In one scRNA-seq sample I encountered the genes: Hbb-bs, Hba-a1 and Hba-a2. These genes appear on top of the list of the highest expressed genes but the 75% percentile of cells have very low counts. Therefore the expression is highly skewed.

Reading from gene cards, they appear to be coding hemoglobin subunits and involved in oxygen transport. I think they may need to be removed before normalization, also because the analysis I'm running does not consider oxygen transport.
Has anyone encountered this situation and what do I do?

You probably have red blood cell contamination.  It's probably better to identify which cells they are.  They'll probably cluster away from the cells you are interested in.


Why ORF13 and ORF14 of Bat Sars coronavirus Rp3 have no corrispondence in Sars2? I'am comparing bats coronavirus with Sars2, i.e. COVID-19. I have one question...

one Bat Sars coronaviurs identify in GenBank as DQ071615.1 has the ORF13(28122..28415) and ORF14(28572..28784) proteins
Sars2 identify in GenBank as MN908947_3, has range 28284..28577
not labeled as ORF13 and the same 28734..28955 is not labeled as ORF14?

The ORF13(28122..28415) of Bat coronavirus and  PSEUDO_ORF13(28284..28577) of sars2 share 72% homology
and ORF14(28572..28784) of Bat coronavirus and  PSEUDO_ORF14(28734..28955) of sars2 share 75% homology

This seems one publication for ORF14 and ORF13 in Bat Viruses and relations in Sars2 here
The author say that ORF14 in Sars2 is all contained in N protein...
And ORF13 is know in that above link as ORF9B
and in that written is said ORF9B is implied in immunitary innate suppression...

The virus Rp3 stands for Rhi pea isolated in China in 2004.
The trees is here and groups with the Bulgarian sarbecoviruses and is closely related to the African sarbecoviruses. Thus although it is a Chinese sarbecovirus is not with the "Chinease originating" clade as is SARS and Covid-19.
Thats why the more mysterious ORFs ORF13 and 14 are missing, because over notable genetic distances these genes along with ORF10 are present/absent. This is a common feature of the betacoronaviruses.
In summmary this virus is notably distance from the major human pathogenic SARS/covid-19 viruses. Thus the absence of these auxilary genes isn't really surprising.
Its the very last taxa on the SARS-related viruses clade (monophyly) and is in a different group to SARS, which is the Beta-1 clade on this most excellent paper.

To answer the comment. They are outgroups, basically ancestral to COVID-19 and SARS-1 clade. So they lost these genes whilst the ingroup including COVID-109 and SARS-1 retained ORF13, 14 ... but I think lost ORF10 from memory


Voom function from limma package and Normalization on counts data I know that Voom function from limma package from Bioconductor converts raw counts into log-CPM values and then Normalization is applied on that, with normalize.method argument.
I would like to know clearly

how this normalization is done (statistically) ?

normalization across samples and normalization across genes. But don't the know how it is done statistically.

I'd like to know answer for this myself. It seems the normalize.method argument is set to none by default. You can edit by looking at the methods in normalizeBetweenArrays.
But how do I know if and when I have to use normalize.method in my analysis?
The output from voom is log2CPM of count values (adding 0.5 to them).


Is re-normalization of RNAseq data recommended for analysis of gene subsets? I downloaded an RNAseq dataset from TCGA database in 3 formats: 1) HTSeq counts; 2) FPKM; 3) FPQM-upper quartile normalized.

The complete dataset contains ~60,000 genes. All of my analysis will focus on a subset of ~2500 genes.

What confused me is the following statement on the TCGA documentation page:
"RNA-Seq expression level read counts produced by HT-Seq are normalized using two similar methods: FPKM and FPKM-UQ. Normalized values should be used only within the context of the entire gene set. Users are encouraged to normalize raw read count values if a subset of genes is investigated."

If I normalize the 2500 gene counts in isolation of the rest, I am guessing the normalization factor will be more unstable between samples than using the complete gene set. It means that if one patient has higher expression of the subset than another, normalization within subset will obscure that difference.

However, the normalization factor (either total mapped reads or upper quartile) becomes more stable the larger the number of genes considered. Less variation will be seen on a global scale because of biological differences, whilst technical differences will be normalized.

In short, do you normalize the gene counts de novo if you analyse only a subset of genes? Why?

FPKM, RPKM are ways of normalising for both sequencing depth and gene length. You first calculate a scaling factor (SF = #read/1Mil),you then divide the genes counts for the SF(RPM = # counts/SF), and finally you divide by the length of the gene (RPM/length(gene)). Each sample has a different SF meaning that you should not compare values from different samples.

To fix this, you have to normalise across samples, by upper quartile normalise using a subset of genes. In this case, the SF is calculated based on the upper quartile (70%) of this shared and immutable subset, and it is therefore the same for all the samples.

If you subset your genes, and you are using the normalised values for downstream analyses, then you need to re-scale your values.

I personally try to avoid FPKM,RPKM,TPM (if I can) and always start from raw counts.

I would suggest to have a look at the documentation of packages like DESeq2 and edgeR for detailed and better explanations.


You could also take a look at Infernal.

I think your best bet is to use RNA-focused programs; I'm not aware of a DNA-specific one, which makes some sense because DNA is typically double-stranded, so the secondary structure community is primarily focused on RNA. I would predict that secondary structure potential of single-stranded DNA would be similar to that of RNA.


WGCNA module preservation I am completing the online tutorials for module preservation of human brain modules in chimpanzee brains and vice versa (III). If anyone is aware of this tutorial or practiced before wondering if can help with quick question please.
I completed sections 1a and 2a to upload the expression data and then proceeded to Appendix A to complete Network Analysis of expression data. Within the code; I get error: Object ConnectivityHuman or Error: Object ConnectivityChimp not found. I have searched for any Sections of tutorial missed but unable to find these? Anyone any advice on tutorial gratefully appreciated thank you

ConnectivityHuman = ConnectivityHuman/max(ConnectivityHuman)
    ConnectivityChimp = ConnectivityChimp/max(ConnectivityChimp)

Thanks for pointing this out. Looks like we forgot to include two important lines. Something along the lines of

ConnectivityHuman = colSums(adjHuman)-1
ConnectivityChimp = colSums(adjChimp)-1

This should do it; I'll update the tutorial.


Align the HIV-1 Gag protein with the Gag protein of Visna virus Align the HIV-1 Gag protein with the Gag protein of Visna virus
Background Visna virus is a lentivirus causing encephalitis in sheep
The problem is align the follow two protein Gag of viruses HIV-1 (GenBank KU749412.1) and Gag Visna virus (GenBank L06906.1) that are:
"MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTTNAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ"
"MAKQGSKEKKGYPELKEVIKATCKIRVGPGKETLTEGNCLWALKTIDFIFEDLKTEPWTITKMYTVWDRLKGLTPEETSKREFASLQATLACIMCSQMGMKPETVQAAKGIISMKEGLHENKEAKGEKVEQLYPNLEKHREVYPIVNLQAGGRSWKAVESVVFQQLQTVAMQHGLVSEDFERQLAYYATTWTSKDILEVLAMMPGNRAQKELIQGKLNEEAERWVRQNPPGPNVLTVDQIMGVGQTNQQASQANMDQARQICLQWVITALRSVRHMSHRPGNPMLVKQKNTESYEDFIARLLEAIDAEPVTDPIKTYLKVTLSYTNASTDCQKQMDRTLGTRVQQATVEEKMQACRDVGSEGFKMQLLAQALRPQGKAGQKGVNQKCYNCGKPGHLARQCRQGIICHHCGKRGHMQKDCRQKKQQGNNRRGPRVVPSAPPML"

I propose that the better alignment method is the one that maximize the function
f(string1, string2)=(number of letter equal between strign1 and string2)/(common length of strings).

With that definition I found f(s1,s2)=0.25760 (and 0.2949 between POL protein) for the two strings above.

What score might others find?
Is that definition one of biological significative?
Is there is a better one, what definition for align proteins?
Are Visna and HIV viruses relatives?

Both are lentiviruses and env is a fast evolving gene. What you are looking at here is the location of indels and this will vary between alignment algorithms.
There will always be an element of error in attempting to align over this genetic distance. It depends what you end goal is, if it is to maximise homology than okay, however optical enhancement of alignments are common if your downstream analysis is phylogenetics. If it is protein structure homology style modelling, where you have the protein structure of HIV env experimentally determined mulitple times you want to use it to help estimate that of the cat lentivirus then optimising a scoring metric is okay.
What you have to understand is there is no 'right' method, in phylogenetics this doesn't matter as much because indels are excluded in the subsequent calculation, so it is less sensative to indel positioning errors and optical enhancement is used. You can also simply exclude regions that appear 'garbled'. For protein homology modeling this is a big problem, because the indels will make drastic changes to the structure prediction. So the only approach here in my opinion is to produce a separate structural prediction for every variation in indel alignment you have obtained. In the latter context getting the alignment right is everything and I suspect thats the output here, because its a classic surface antigen to understand, which will readily crystalise.


Minimap2 -ax map-pb doesn't output tlen field I have used minimap2 to map some pacbio reads to a reference genome. I would like to know the "insert size" (true length of sequence) relative to the reference. More specifically, I want to know the first base position in the reference where the read maps and the last base position in the reference where the mapping ends.

As I understand it I need the TLEN field from bam files (sam format field 9). However, all the alignment files from minimap2 have 0 in the TLEN field.

Does anyone know how to get pacbio alignments with the TLEN field? I am also happy using bwa or another aligner.

The command I used was minimap2 -ax map-pb {fasta} {reads} | samtools view -h -b - | samtools sort - > myfile.bam.

The TLEN for single-end reads is defined as 0 in the SAM specification, which is why this happens. You can still get the end position of the alignment by parsing the CIGAR string. An example using the pysam package in python is below:

import pysam

bam = pysam.AlignmentFile("myfile.bam")
for b in bam.fetch():
    print("This alignment spans {}:{}-{}".format(b.reference_name, b.reference_start+1, b.reference_end))

The reference_end value is automatically calculated from the CIGAR string.


How to calculate AUC in coverage graph Is there a way to calculate the area under a curve in a coverage graph?

Thank you in advance.

EDIT:

I'd like to calculate the AUC per strand per gene in a coverage graph. I have two wig files (sense and antisense) generated from a single bam (RNA-seq) file (using the strand-coverage software).

Yes, just add up the coverage at each location; that's integration in a nutshell.

Or, alternatively / equivalently, count the total number of reads within the region.


How to edit the headers of multiple fasta files from multiple folders? (recursively) This is an expansion of my previous question, How to edit the headers of multiple fasta files from multiple folders?

My directories are organized as follows: one main directory, in which I have multiple project directories. Inside each project directory, there are multiple metagenome directories. Inside each metagenome directory, there is always a single .fna (fasta) file.

Names: project directories start with 'mgp', metagenome directories start with 'mgm' and end with '.3'. They end in numbers, but the numbers are not related. (eg. project directory 'mgp83581' contains metagenome directory 'mgm4729322.3').

What I want is a loop that would go through all these fasta files, and add to the headers codes corresponding to project and metagenome directory. For example, if the fasta file in project directory mgp83581, metagenome subdirectory mgm4729322.3 had the following headers...

>seq1

>seq2

...then I would want to change it to:

>83581_322_seq1

>83581_322_seq2

Basically, I want to add the number code associated with the project, and the last three digits of the 7-digit code associated with the metagenome (apart from the '.3' part of the file name).

Below is the answer that I got for my previous question, where I didn't have project directories, and where the 3-digit code for the metagenomes stayed in a manageable range. (thanks terdon!)

for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done

Using ubuntu. I would appreciate your help!

You can edit your for loop so as to include the paths to your files of interest.

The command below will give you the relative paths (relative to where you run the command, I think it is better to run this within your main directory) of your files of interest:

ls mpg*/mgm*.3/*.fna

You can use the output of this command (your files of interest) to iterate over via your for loop. Your code will look like:

for i in {...}; do
    sed -i.bak "s/^>/>${i}_/" i
done


Merging regions according to their identifier I have a huge file (20 GB) which has a range of genomic locations, and for each location there is an identifier(4th column), which is sometimes the same.

file1.txt

chr1  10  20  ABC
chr1  13  20  ABC
chr1  14  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  21  JJJ
chr2  21  31  JJJ
chr2  23  27  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  101 MMM
chr3  101 105 MMM

I am trying to merge the ranges of column 2 and 3, if they are in the same chromosome and having the same identifier ( having the 1st and 4th columns the same).

I have tried to use bedtools, as shown in this post, which looks very similar but not the same with mine, as they are looking for different outputs ( range vs merge )

So the answers in that post was with using groupby in bedtools2,source code from github.

I have applied the same as the following:

./groupBy -i ../file1.txt -g 1,4 -c 2,3 -o min,max | awk -v OFS='\t' '{print $1, $3, $4, $2}'

chr1    10      37      ABC
chr2    15      56      JJJ
chr2    25      26      MMM
chr3    32      105     MMM

But it only connects the ranges of identifiers.. It groups by, and does not merge.

And, if I use bedtools v2.26 merge function :

>sort -k1,1 -k2,2n file1.txt > file2.txt
>cat file2.txt
chr1    10      20      ABC
chr1    13      20      ABC
chr1    14      21      ABC
chr1    22      27      ABC
chr1    29      37      ABC
chr2    15      21      JJJ
chr2    21      31      JJJ
chr2    23      27      JJJ
chr2    25      26      MMM
chr2    35      56      JJJ
chr3    32      42      MMM
chr3    45      76      MMM
chr3    88      101     MMM
chr3    101     105     MMM

bedtools merge -i file2.txt
chr1    10      21
chr1    22      27
chr1    29      37
chr2    15      31
chr2    35      56
chr3    32      42
chr3    45      76
chr3    88      105

Which, merges without seeing the identifiers as it also gives me this output.

What I want is to merge the ranges if they overlap, if having the same 1st and 4th column as:

expected_output.txt

chr1  10  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  31  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  105 MMM

I only have access to 15GB of RAM, so I can't store the entire file in memory.

If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g.

awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \
| mergeBed \
| sed 's/_/\t/' \
| awk -v OFS='\t' '{print $1, $3, $4, $2}'

(Assuming "_" is not part of chromosome or identifier name, adjust accordingly if so).


What is the typical host-to-bug DNA ratio found in nanopore sequencing without amplification? I'm interested in sequencing a sputum sample using an ONT MinION without performing any type of DNA amplification.
Has anyone found a good reference (or anecdotal evidence) for a range of human to pathogen DNA ratios typically found in sputum samples?
Does this depend on whether the patient is infected by a virus or bacteria, and the bacterial species?
Can I expect to find 10x as much human DNA as pathogen DNA, or 100x-1000x? I would really just like a ballpark estimate.

My recollection of a conversation I had with Justin O'Grady is that it's about 99% host DNA for unfiltered samples. There are a few ways to reduce that portion prior to sequencing; he has written a paper that mentions a saponin-based host depletion method that is quite effective:
https://doi.org/10.1038/s41587-019-0156-5


How to link GDC ids to CCLE cell line names? Hello bioinformaticians,
I've recently downloaded few RNAseq data files from Genomic Data Commons Data (GDC) portal. These files belong to Broad institutes CCLE project.

Now the problem is that GDC portal do not provide CCLE cell line names:

So my question is how could I link the GDC ids (prefixes of File Names) to CCLE cell line names?

I would imagine because it is an important step in the analysis, there should be some way to do that.

The filenames of those bams contain the cell line names. They are constructed:

SAMPLEID.CELL_LINE_NAME.REPLICATE.bam

Thus the first 5 cell lines in your screen capture are SK-N-F1, TE-9, CJM, HCI-H1915 and CAL-33.

All the information that is need is also present in the experimental metadata xml file attached to each sample record. This is accessible by downloading it from the record on the website, or by querying the GDC metadata API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api).


Question on nanopore sequencing data process pipeline (cDNA-PCR) I recently started doing the analysis on nanopore sequencing data. As I was searching for some help on pre-processing of the data, I found your nice setup pipeline created here: https://www.protocols.io/view/stranded-mapping-from-long-reads-n8cdhsw.

As I was trying to follow the instruction, I ran into a question on part where it uses lastal to divide the input file as strand specific output.

Part that uses the following script below:

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \   awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \   ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \   awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

In here, I was wondering what the ONT_barcodes_adapters.fa file actually looks like. I wanted to find the information in other places and could not find it elsewhere.

If you could help me with this, I would really appreciate it!

I recommend that you look at the most recent version of the protocol, which includes example adapter and barcode files:

https://www.protocols.io/view/stranded-mapping-from-long-reads-syheft6?step=2

[check in the protocol to see if there's a more recent version]

These files contains the public sequences that ONT have released for their sequencing kits (as briefly mentioned in the "Before Starting" section of the protocol). The sequences are available on the ONT community, but I have converted them into FASTA files for use in my protocol.

The adapter and barcode files were combined into a single FASTA file in the previous version of my protocol. For the PCR cDNA barcoding, it looked something like this:

>ONT_VNP
ACTTGCCTGTCGCTCTATCTTCTTTTT
>ONT_SSP
TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG
>BC01
AAGAAAGTTGTCGGTGTCTTTGTG
>BC02
TCGATTCCGTTTGTAGTCGTCTGT
>BC03
GAGTCTTGTGTCCCAGTTACCAGG
>BC04
TTCGGATTCTATCGTGTTTCCCTA
>BC05
CTTGTCCAGGGTTTGTGTAACCTT
>BC06
TTCTCGCAAAGGCAGAAAGTAGTC
>BC07
GTGTTACCGTGGGAATGAATCCTT
>BC08
TTCAGGGAACAAACCAAGTTACGT
>BC09
AACTAGGCACAGCGAGTCTTGGTT
>BC10
AAGCGTTGAAACCTTTGTCCTCTC
>BC11
GTTTCATCTATCGGAGGGAATGGA
>BC12
CAGGTAGAAAGAAGCAGAATCGGA
>RB12A
GTTGAGTTACAAAGCACCGATCAG

I separate out adapters and barcodes in the newer versions of the protocol because I was getting too many incorrect hits.


Mapping Reads to Known Gene Paralogs with Long Read Technology I have some sequencing data from a captured region that is a known paralog edited. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:

minimap2 -ax map-pb GRCh37.fa fastqfile.fastq > samfile.sam

I captured a 10kb region and what I am seeing in IGV is puzzling. The region directly before and directly after the paralogues are covered at about 1000X, but the genes themselves are not covered at all:

The two genes in the middle are the paralogues. I left them off to protect data. This data is also only the sorted data, the duplicate marked data has lower coverage (600X).

Question 1

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Question 2

Is this process normally accomplished by minimap2, or are other aligners more effective for paralogous mapping?

Question 3

The interesting thing is, while my average read length was 6kb, I am seeing the reads here are about 4kb. I have another region with very high coverage where the reads spanned both paralogues and mapping was successful.

Question: Why am I only seeing reads that are about 4kb long when the average read length for this dataset is 6kb?

Any help on this would be great!
Thanks and I hope I have been clear on my question,

Dennis

Question 1

This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.

samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam

This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the bam file.  Using picard tools is the fastest way I have seen to filter a bam file given a list of reads.

#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$filtint" ${BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq > ${TEMP}; \
java -jar ${PICARD_JAR} FilterSamReads \
  I=${BAMFILE} \
  O=${OUTFILE} \
  READ_LIST_FILE=${TEMP} \
  FILTER=includeReadList
rm ${TEMP}

This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.

Question 2

Filtering the bam file should be pretty standard, as addressed above in my response to Q1. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.

The fact that there are paralogs does not matter to the aligners. The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. I do not know of any long-read paralog-aware aligners.

As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.

Question 3

Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.

One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.


What are common ways to calculate gene length for TPM calculations? I am aware of this similar question. But the accepted answer there answers how to calculate TPM given a mapping from gene name to gene length.
My question is, given an annotation file of genes (for example from Ensembl here), how do you calculate gene length for the purpose of calculating TPM on a gene count matrix coming from RNA seq data. That is, I understand the actualy calculation of TPM given gene lengths, but I do not understand what the gene lengths are.
From what I understand, to calculate gene lengths from annotation file for RNA seq data, you only consider the exons of each gene, but even so, many exons overlap. What is the standard method of dealing with this?
For example, I have read that using the length of the union of all exons of the gene is acceptable, but from what I understand TPM is meant to represent the "expected number of reads coming from gene X when randomly selecting 1e6 transcripts". So taking the union of all exons will be an estimate of this but it might under-estimate this.
EDIT: From what I can tell this question is not a duplicate of RNASeq: Normalization, stabilization, gene length and rlog
That question does discuss how to perform normalization that takes gene length into account, but it does not ask "what is gene length? How does one determine the length of each gene?" Possible answers could be like what I already mentioned: length of the union of all exons in the gene, or another possible answer is one from the linked question "length of the longest transcript isoform" but I feel that it is not really discussed in that answer/question.

I think the most common way is to use something like Kallisto which figures out what proportion of reads are likely assigned to each transcript variant, and calculated TPM  based on that.