text:
How to plot character state changes from a presence-absence matrix on to a phylogeny
I wish to assign character state changes from a presence-absence matrix to a phylogeny. Thus I want to identify the most parsimonious hypothesis for which node and branch a given "mutation" or "change" has occurred.
I have tried assigning each character to its leaf node, and then if the leaf node's sister has the same character I reassign the character to the parent node (and work back until all nodes are assigned). I am using a dummy dataset to try to achieve this:
Matrix
>Dme_001
1110000000000111
>Dme_002
1110000000000011
>Cfa_001
0110000000000011
>Mms_001
0110000000000011
>Hsa_001
0110000000000010
>Ptr_002
0110000000000011
>Mmu_002
0110000000000011
>Hsa_002
0110000000000011
>Ptr_001
0110000000000011
>Mmu_001
0110000000000011

Phylogeny
((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));

I assign internal nodes using ete3, so my output should be:
BranchID    CharacterState    Change
Node_1:    0    0->1
Hsa_001:    15    1->0

As my code assigns character states based on their sisters if a loss is encountered it messes up the output so that:
BranchID    CharacterState    Change
Node_1:   0     0->1
Node_3    15    0->1
Node_5    15    0->1
Node_8    15    0->1

I'm coding in python 2.7 and welcome assistance.
My code:
from ete3 import PhyloTree
from collections import Counter
import itertools

PAM = open('PAM','r')

gene_tree = '((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));'

NodeIDs = []

tree = PhyloTree(gene_tree)
edge = 0
for node in tree.traverse():
    if not node.is_leaf():
        node.name = "Node_%d" %edge
        edge +=1
        NodeIDs.append(node.name)
    if node.is_leaf():
        NodeIDs.append(node.name)

f = open('PAM','r')

taxa = []
pap = []

for line in f:
    term = line.strip().split('\t')
    taxa.append(term[0])
    p = [p for p in term[1]]
    pap.append(p)

statesD = dict(zip(taxa, pap))

def PlotCharacterStates():

    Plots = []

    events = []

    for key, value in statesD.iteritems():
        count = -1
        for s in value:
            count+=1
            if s == CharacterState:
                a = key, count
                events.append(a)

    Round3_events = []
    while len(events) > 0:
        for rel in Relationships:
            node_store = []
            sis_store = []
            for event in events:
                if rel[0] == event[0]:
                    node_store.append(event[1])
                if rel[1] == event[0]:
                    sis_store.append(event[1])
            if (len(node_store) > 0) and (len(sis_store) > 0):
                place = rel, node_store, sis_store
                Round3_events.append(place)

        moved = []
        for placement in Round3_events:
            intercept = (set(placement[1]) & set(placement[2]))
            node_plot = (set(placement[1]) - set(placement[2]))
            sis_plot = (set(placement[2]) - set(placement[1]))
            if len(node_plot) > 0:
                for x in node_plot:
                    y = placement[0][0], x
                    Plots.append(y)
                    moved.append(y)
            if len(sis_plot) > 0:
                for x in sis_plot:
                    y = placement[0][1], x
                    Plots.append(y)
                    moved.append(y)
            if len(intercept) > 0:
                for x in intercept:
                    y = placement[0][2], x
                    y1 = placement[0][0], x
                    y2 = placement[0][1], x
                    moved.append(y1)
                    moved.append(y2)
                    events.append(y)

        for event in events:
            if event[0] == "Node_0":
                Plots.append(event)
                moved.append(event)

        events2 = (set(events) - set(moved))
        events = []
        for event in events2:
            events.append(event)

    pl = set(Plots)
    Plots = []
    for p in pl:
        Plots.append(p)

    print CharacterState, Plots

'''
assign sisters to leaves, internals
'''

e = []
round1b_e = []
round2a_e = []
placements = []
Relationships = []
Rounds = []
for node in tree.traverse():
    sisters = node.get_sisters()
    parent = node.up
    cycle1 = []
    if node.is_leaf():
        for sister in sisters:
            if sister.is_leaf():
                round1a = ["Round1a", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1a)
                e.append(node_names)
                x = node.name, sister.name, parent.name, "leaf-leaf"
                Relationships.append(x)
            if not sister.is_leaf():
                round1b =  ["Round1b", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1b)
                round1b_e.append(node_names)
                x = node.name, sister.name, parent.name, "node-leaf"
                Relationships.append(x)
    elif not node.is_leaf():
        if not node.is_root():
            for sister in sisters:
                if not sister.is_leaf():
                    node_names = node.name, sister.name
                    round2a_e.append(node_names)
                    x = node.name, sister.name, parent.name, "node-node"
                    Relationships.append(x)

x = []
CharacterStates = []
for key, value in statesD.iteritems():
    for value in value:
        x.append(value)

y = sorted(set(x))
for x in y:
    CharacterStates.append(x)

for CharacterState in CharacterStates:
    PlotCharacterStates()

My code is added, thank you for the suggestion. I am using python 2.7, and I am not yet using biopython but happy to use it if there is a solution. I am looking to replicate the -apo function from TnT (Goloboff 2008) if possible (if you are familiar with this)

Please [edit] the question to include this information, it might be helpful for others (I am not familiar with these area, sorry).

To give some formal words for what you are doing, I would say this sounds like you are trying to reconstruct ancestral states.  There has been some work done on this, and a similar algorithm to yours was devised by Fitch in the 1970s: Fitch W: Toward defining the course of evolution: minimum change for a specific tree topology. Syst Zool 1971, 20: 406–416. 10.2307/2412116
The missing part, from your approach, is that once Fitch gets to the root of the tree, he than retraverses back to the leaves, resolving ambiguities along the way (if two sister nodes have different states, how do you decide what the parental state was? Fitch's idea is that once you get to the root, you can come back down and use the parental state to set the child state in ambiguous cases).
There is a nice implementation of Fitch's algorithm in the R phangorn package, which makes the reconstruction a one-liner (given a rooted tree object, and a matrix of leaf node states):
anc.mpr = ancestral.pars(rooted_tree, tree_samples_data_matrix.phyDat, "MPR")

If you don't want to use the R package itself, you could perhaps use it as pseudocode to port into python.

The OP could use the Rpy2 library to connect Python to R, but it is often better to use R directly.
The alternative is to simply use MacClade 4.0, although you will need an old Mac. This is a GUI interface that will allow you to identfy the parsimonious character state for each node in the phylogeny. The API allows you to interactively explore the phylogeny by point and click moving the branches around the tree. The tree will then be assessed in terms of "number character steps" and automatically calculate this for each topology loaded. The output is very colourful so its easy to understand a large data set.


text:
What is the index fastq file (sample_I*.fastq.gz) generated when demultiplexing Illumina paired-end runs?
What is the index fastq file that comes with some Illumina sequencing datasets? (The samplename_I*.fastq.gz file.)
For example, I recently received some 10X Chromium reads for two libraries sequenced on the same lane. This was a 2x150 sequencing run, so there should be two fastq files. The sequencing center demultiplexed the libraries and generated two separate directories - one for each library.

In each directory there are three fastq files:

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

I know that the R1 fastq file contains the forward reads and the R2 fastq file contains the reverse reads, but what is the index file? The first few sequences do not match the 10X indices for this library (ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA).
>bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | head
NTTGGGTA
NGCAATGG
NAGCCCAC
NAGCCCAC
NCATTACT
NCATTACT
NGCAATGG
NAGCCCAC
NTTGGGTA
NAGCCCAC

I don't know 10x Chromium software, but I know their single cell pipeline, and many versions ago, the software was expecting the read indices to be in a separate file along with the reads themselves.  Cellranger has since changed their pipeline so that it doesn't want these anymore, but whoever gave you your files might have seen that it was a 10xGenomics protocol, and decided to give you the indices like that just in case you needed them,

For what it's worth, 10X doesn't require the index files as input (you can just give it demultiplexed samples) and Illumina's software won't produce those index files by default. Also, you can have up to 2 index files for Illumina data (though you typically won't if you're making standard 10X libraries).

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

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

It took me a while to figure out that the "index" is the same thing as the "barcode" that says which sample each sequence is from on a multiplexed run.
If your data is not demultiplexed (single R1.fastq and R2.fastq files contain the information for multiple samples), then this I1.fastq file is what you use to assign each sequence to a sample (ie to "demultiplex"). The spreadsheet listing which barcode/index is which sample is often called a "mapping file."
If your data is already demultiplexed (you have separate R1.fastq and R2.fastq files for each sample), then you don't really need it... however some people use the index sequences in their quality control workflows. For example, this 2016 paper in BMC Genomics by Wright and Vetsigien "Quality filtering of Illumina index reads mitigates sample cross-talk": https://doi.org/10.1186/s12864-016-3217-x

The term index is almost always used to refer to the 8-bp sequence that the Illumina machine reads to identify the library-of-origin of the individual read, while the term barcode could also mean some sequence within R1 or R2 itself that the user added to label provenance - such as being from a specific cell or specific amplicon.

Yes, in the context of a 10xGenomics single cell run, the "barcode" means the cell barcode, which is not at all the same thing as the sample index.

I received 3 files that exactly match the file names listed by the original poster, except that my data was clearly not demultiplexed yet. My index file referred to what my mapping file listed as the "SampleBarcode".  It's certainly possible that formats are different for single-cell data. For perspective, my data is from a multiplexed 16S rRNA gene amplicon miSeq run. My barcodes identified which sample each sequence was from, and they were a 12 bp sequence attached to the forward primer (so it was not dual-indexed). This is the protocol used by the Earth Microbiome Project.

I looked at 10x data a couple of years ago, and at that time, their protocol would use four separate indexes, with balanced index sequences, as a way of introducing diversity into the library.  These four indexes were separate from the individual oligos used to mark single cells.  For an individual sample, you would need to combine the reads from all four separate indexes, which is how it looks like in your data.  There is a brief explanation here: https://support.10xgenomics.com/single-cell-gene-expression/index/doc/specifications-sample-index-sets-for-single-cell-3

source: Stack Exchange

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

sort: votes


text:
Why are my Chi-squared test results different from those in a published table?
I recently read the paper “A novel long non-coding RNA linc-ZNF469-3 promotes lung metastasis through miR-574-5p-ZEB1 axis in triple negative breast cancer”. In this I see Table1 showing correlation of Linc-ZNF469-3 with different features in TNBC and non-TNBC patients. p-values are calculated based on Pearson's Chi-sq test.

For non-TNBC patients: I calculated the p-value with the Chi-sq test just to have a look if everything is right.
Recurrence:
data <- data.frame(x= c(40,62), y=c(56,65))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 1.1269, df = 1, p-value = 0.2884

Brain:
data <- data.frame(x= c(6,99), y=c(6,119))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.096461, df = 1, p-value = 0.7561

You can see that the Recurrence p-value doesn't match with what I see in the table. Others, like Brain, Bone, Liver, Lung match with the p-values I calculated.
And for TNBC patients, not a single p-value matches with my calculations. I don't even get Lung and Bone as having significant p-values with my calculations.
I did it for Lung as shown here:
data <- data.frame(x= c(25,143), y=c(9,35))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.8044, df = 1, p-value = 0.3698

But in the paper they mentioned:

We found that high linc-ZNF469-3 expression was significantly correlated to lung metastasis in TNBC patients (Table 1).

Can anyone tell me what could be the wrong here? Is it my calculation wrong? or in the paper?

Ouch. This probably warrants a correction or even retraction.

Exactly. I have another question. Lets say Lung p-value in TNBC patients is significant. High linc-ZNF469-3 has 9 and low linc-ZNF469-3 has 25. How could one tell that high linc-ZNF469-3 expression was significantly correlated to lung metastasis in TNBC patients?

We found that high linc-ZNF469-3 expression was ignificantly correlated to lung metastasis in TNBC patients (Table 1). This conclusion  can not be reached. Metastasis across tnbc and non-tnbc require a different analysis e.g. t-test.

Your calculations seem right, perhaps there was an error on their side.
I also looked at the number of samples reported, but they use the same amount of samples in each case. Because they are testing whether Linc-ZNF469-3 is different in several tissues, they should have used a multiple test correction for the p-value (although reporting the original values is usually also ok).
To answer your question in the comments. You can't correlate a categorical variable with a continuous  variable. Or they did another test or they are not explaining clearly what they did.

Thanks for the reply. For my question...I also had a look into some other papers. For eg: Check this table [https://www.nature.com/articles/s41598-017-04736-z/tables/1] They say in the paper like this; High RAD51-AS1 expression correlated with advanced FIGO stage, histological subtype, high histological grade, number of tumors and distant metastasis. I'm bit confused. Is it based on their percentage they say RAD51-AS1 expression is correlated with those variables?

@user3351523 They could show a plot of like FIGO stage and the expression of the variable, and see a linear trend towards higher expression on certain stages, but without comparing the mean or the percentage/proportion (and not a correlation) that could be noise or illusions.

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


text:
How do I rewrite a read group using pysam?
I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups?
import argparse
import pysam

def rewrite_read_groups(input_filename, output_filename, sample):
    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        # replace header fields
        output_header = bam_input.header
        read_groups = output_header['RG']
        print(read_groups)
        for read_group in read_groups:
            read_group['SM'] = sample
            print(read_group)
        output_header.to_dict()
        print(output_header)
    finally:
        bam_input.close()

    print(output_header)

    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)

        # copy alignments
        for read in bam_input.fetch(until_eof=True):
            bam_output.write(read)

    finally:
        bam_input.close()
        bam_output.close()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Rewrite bam read groups")
    parser.add_argument('-s', "--sample", help="Set sample for all read groups")

    parser.add_argument('-i', "--input_file", help="Input bam", required=True)
    parser.add_argument('-o', "--output_file", help="Output bam", required=True)

    args = parser.parse_args()

    rewrite_read_groups(args.input_file, args.output_file, args.sample)

If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups.
Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not to be stable across PySAM versions.

I want to assign all read groups to the same sample, where read groups are from different flowcells, but not rewrite all alignments to the same read group, so the Picard AddOrReplaceReadGroups is not the tool I am looking for. I have actually now written something using htsjdk in java to perform the operation I want, but I will leave the question for pysam.


text:
Why are there more barcodes than GEMs in 10X chromium data?
Question: Why are there more barcodes than GEMs in 10X chromium data?
Introduction
I have several 10X genomics chromium libraries made with the de novo assembly/genome protocol. The 10X chromium device is advertised to produce approximately 1,000,000 GEMs (gel-emulsion droplets). In theory each GEM should have a single barcode that is used to track the provenance of each read pair. If this is true, a single library should have approximately one million barcodes, whether or not it was overloaded or underloaded.
Results
The libraries actually have more than one million barcodes, confirmed by longranger basic. Library 1 has 2,369,949 barcodes. Library 2 has 2,311,405 barcodes. Both libraries have approximately 150,000,000 read pairs.
Why this is a problem
I'm using the 10X data for some non-standard analyses and want a high number of LinkedReads per molecule (LPM). Having multiple barcodes in a GEM will reduce the number of LPMs per barcode proportional to how many barcodes were present in the GEM.
Possible explanations

If there are more than approximately one million barcodes in a library then this leaves several possibilities:

Each GEM truly has a single barcode. Therefore the Chromium device actually produced significantly more GEMs than advertised.
Each GEM has a variable number of barcodes. Therefore the number of barcodes is independent of how many GEMs are created.
Some combination of the above two options. Some variable number of GEMs and a variable number of GEMs per barcode.
Additional barcodes are created by sequencing errors and misclassification by longranger basic. This would probably manifest in many barcodes that only have one or two linked reads total.

Deeper analysis
I don't know how to test hypotheses 1, 2, & 3 bioinformatically. Hypothesis 4 is testable though by making a histogram of read counts per barcode. To do so, we will use the output of longranger basic for Lib2.
Longranger basic outputs an interleaved fastq.gz file that contains the whitelisted barcode in the header. For example (pretending the reads are shorter than 150bp):
@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
AATTACCCAAAATTAATATTTGTATTAATTTTCATATTAATAAC
+
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
CACAATATCAAAAACAAATCAAAAATTAAGAAAATATAACTAAG
+
AAF<7FF-JJFJJFJJJJJJAJFJJJFJJJJJFJJJJJJJJJJJ

The BX:Z:AAACACCAGACAATAC-1 tells us that the index is whitelisted (BX:Z), the index sequence (AAACACCAGACAATAC), and the library index for the analysis (-1).
So, we can count the number of occurrences of each barcode with this code:
bioawk -cfastx '{print($comment)}' barcoded.fastq.gz | sort | uniq -c > bc_counts.txt

The first column of bc_counts.txt tells us how many individual reads had that barcode. We divide this number by two to see how many read pairs had that barcode and make a histogram.
awk '{print($1)/2}' barcode_counts.txt | sort | uniq -c | sort -k2 -n | awk '{printf("%d\t%d\n", $2, $1)}' > pairs_per_bc.txt

The first column of this file has the number of read pairs in the barcode, and the second column has how many barcodes have that number of read pairs. We see that there are many barcodes that only have one or two read pairs:
1       593583
2       207192
3       79409
4       34174
5       17353
6       10662
7       7997

It looks like every barcode with fewer than 10 read pairs may be noise. How many barcodes contain at least 10 reads?
>awk '{if ($1 >= 10){sum += $2} } END{print(sum)}' pairs_per_bc.txt
1347686

Around one million.
Conclusion
This analysis doesn't give conclusive evidence for the four hypotheses above, but it seems like there is a noise-generating process in 10X data that causes there to be many barcodes with only a few readpairs.

This seems a nice study, but what is the question? And what type of answers are you looking for?

I mostly want to know why there are more barcodes than 10X advertises there are GEMs in a single chromium library. The kind of answer I'm looking for is for someone to clarify if there are multiple barcodes per GEM, more GEMs than advertised, or a mixture of those two things. Also, it would be nice to know why there is such variation in the number of linked reads per molecule.

Your initial guess is almost certainly correct. I don't know about the linked read libraries, but in the 10X single cell sequencing protocol, separating real barcodes from noise barcodes is an important and sensitive step in the analysis pipeline.
The current analysis pipeline is to look at a plot, like the one you produced and find the "knees". The plot you produced has two knees, on at around 10 and one at around 200. I'd guess that the ones below 10 are what we call "error barcodes" generated by sequencing is PCR errors from real barcodes, and the ones between 10 and 200 are what we'd call "background" barcodes - in scRNAseq they represent things like GEMs that didn't encapsulate a cell, but did encapsulate some non-cellular DNA or RNA.
As there are more than 1,000,000 possible barcodes, not all of them are used.  I believe 10X knows which sequences are used and which are not (ie the barcodes are not random), so I'm surprised that these sequences make it out of the other end of the analysis.


text:
finding unique headers in a fasta file using linux command line
I tried to use the following command
uniq -u reference.fasta >> reference_uniq.fasta

I'd like a count of the unique headers.

does it matter what the output order of the unique headers is?

If you just want the number of unique headers, you can do this:
grep '>' reference.fasta | sort | uniq | wc -l

If you want a list of the unique headers, you can do this:
grep '>' reference.fasta | sort | uniq

If you want a histogram of how many times each header occurs, you can do this:
grep '>' reference.fasta | sort | uniq -c | awk '{printf("%s\t%s\n", $1, $2)}'

The uniq command expects sorted input. Interestingly, the sort command actually has a "unique" option, -u, which means uniq is not strictly needed. For the fastest processing, you can look for the '>' character at the start of lines with grep:
grep '^>' reference.fasta | sort -u > reference_headers_unique.fasta

For returning the number of unique lines, pipe through wc -l:
grep '^>' reference.fasta | sort -u | wc -l

For more information about regular expressions, see here.

This seems like the best solution given, just note it assumes that the desired output order for the unique headers is the order provided by sort (which may not be stable depending on the locale)

If, for some reason, you need to keep the original order of the input file, so you only want the first occurrence of each header, you can do:
awk '/^>/ && !a[$0]++' reference.fasta > headers

Or, to get the number only:
awk '/^>/ && !a[$0]++{k++}; END{print k}' reference.fasta

However, these only make sense if you need to keep the order of the original file. If not, use Gringer's grep approach which will be much faster (reference.fasta is hg19):
$ time awk '/^>/ && !a[$0]++{k++}; END{print k}' reference.fasta
93

real    0m34.510s
user    0m31.861s
sys     0m2.634s

$ time grep '^>' reference.fasta | sort -u | wc -l
93

real    0m16.597s
user    0m0.820s
sys     0m3.725s

+1 for order capturing approach and benchmark, I guess ripgrep might be even faster


text:
convert tRNA coordinates Result to Fasta format
I believe many of the experts here have gone through this issue,
I have run tRNAscan-SE on my genome to predict all the tRNAs, and have successfully gotten all the tRNA coordinates like start, end and type all  arranged in a table.
My question is how can I convert this result to fasta format as long I have all the details about the tRNA ready?
Is there any tool that I can use to get the fasta sequences of the regions described in my table from the original genome?

Sequence                tRNA    Bounds  tRNA    Anti    Intron Bounds   Inf
Name            tRNA #  Begin   End     Type    Codon   Begin   End     Score   Note
--------        ------  -----   ------  ----    -----   -----   ----    ------  ------
NW_011590950.1  1       36704   36632   Ala     AGC     0       0       40.8    pseudo

Could you provide the table you are referring? Having some sample data will help us to answer your questions. Do you have the reference genome? Many thanks! BTW, what have you tried? Perhaps manually finding the sequence or something else?

I have uploaded the table, unfortunately, the data I'm working on doesn't have a reference.

finding these results just by running tRNAscan on my genome, the table above is my Results, in fact, I need a fasta format has all the tRNA sequence to implemented in another tool

You provided a picture, could you provide text that could be used to copy and paste to make it easier for people to use the information?

I did upload a copy from the data but, when past here doesn't arrange accordingly

I could edit and arrange the data, but consider pasting it and format it using code formatting (that will keep the structure of the data)

What file did you run tRNAscan on? Do you have a single, multifasta file with the sequences NW_011590950.1, NW_011590951.1, NW_011590952.1 etc?

tRNAscan-SE, when run locally, should also provide you with a FASTA output file, if I remember correctly.

this is my command, im runing trnascan locally --> tRNAscan-SE -E ../GCDROM.fna -o MytRNA -f MytRna-secondary -mMytRna-statistics                          but the run end up with an errore in phase 3

Status: Phase I: Searching for tRNAs with HMM-enabled Infernal
Status: Phase II: Infernal verification of candidate tRNAs detected with first-pass scan
FATAL:  Unable to open output file MytRNA

You can select “Output BED format” on the tRNAscan-SE server:

Given the output BED format and your input reference in FASTA format, you can easily generate a FASTA file of tRNA sequences via bedtools getfasta:
bedtools getfasta -name -s -fo ‹output.fasta› -fi ‹reference.fasta› -bed ‹tRNAscan-SE-output.bed›

thank you for the respond, but I'm running it locally, and the output options are Secondary, statistics, and shown table. but currently, I have submitted the job to their server, let's wait and see :)

You can do this using the tools from the exonerate package. I assume you have all the sequences shown in your table in a single, multifasta file. If so, first create a single file per sequence using fastaexplode:
fastaexplode reference.fa

This will create one file for each sequence in reference.fa. Now, you can read through the results table, and extract the sequences:
tmpFile=$(mktemp)
awk -F"\t" '{
                if(NR>3){ printf "%s %s ", $1, $2;
                    if($4>$3){
                        print $3, $4-$3
                    }
                    else{
                        print $4,$3-$4,"-"
                    }
                }
            }' file | while read name num start len rev; do
                fastasubseq "$name".fa $start $len > $num.fa
                if [[ ! -z $rev ]]; then
                    fastarevcomp $num.fa > $tmpFile &&
                    mv $tmpFile $num.fa
                fi
                done

I would suggest you go with Konrad's approach, however, that looks much cleaner and simpler.


text:
Tool to remove a PCR contamination in NGS data
BACKGROUND
I work on NGS data (illumina paired ends reads) coming from a full extract of RNA (metagenomic). We are interested in the viral fraction of this extract.
I observed a contamination with a PCR amplicon. This contamination has been identified as coming from a PCR result that we used to make near the library preparation. So I know exactly the primers used in the PCR.
Here you can see a tablet visualization illustrating the contamination (2 contaminations). I obtained that image with a mapping on a known reference of my virus.
PROBLEM
The problem I have is that these PCR products are exactly the type of sequence I want from my extract (identifying part of the virus). So I have some "false positive" and some misidentification of viruses present in the extract.
I have thought about removing all reads containing primer sequences but it leads to a gap in my data and I can't identify the viruses in my extract anymore.
I am thinking that all of the read pairs that begin before or end after my primers are "true-positives" and I would really like to keep them (even having the primer sequence in them).

Do you know any tool that can help me in that situation? i.e. removing reads containing a sequence but only if it begins with it.
If the tool can also handle the degeneration of my primer it would be perfect.
An example of primer I used is
CNTGGGAGGGCGATCGCAA

Complementary information (asked in comments)
I have around 30 million read pairs (60 millions reads in total). Contamination reads are between 200.000 to a 1 millions reads (depending on sample). "true" reads (not from the contamination) represent around 0.5% of the contamination.

Welcome to the site Untitpoi. Could you please clarify how many read do you remove when you "remove all reads containing primer sequences"? Also how many reads you have in total? Could you quantify the contamination (i.e 50% of your reads are contaminated)? Also how did you try to remove these sequence did you use a bash script or something like that?

Have you searched this site yet?  There are other similar questions such as https://bioinformatics.stackexchange.com/questions/5163/extracting-sequences-from-fasta-beginning-with-common-5-end.   Also what tools have you tried?

Can you post an example of what you reasonably know to be a PCR contaminant and a reasonably similar read that you believe not to be but gets removed using the methods you've tried? That would help us fully grasp exactly what you're dealing with.

@DevonRyan I added some graph to explain my problem, I hope this is clearer now. I didn't try anything yet but I will try a cut-adapt method today and the bioawk method proposed by conchoecia.

I just saw your post with the exact sequence of the forward and reverse primers. Is it true that you expect the beginning of one read to start with the forward primer, and the beginning of its paired read to start with the reverse primer? If so, I should change my answer to reflect that and to better filter your data.

This bash script will get you all of the reads that do not start with your primer sequence. It also handles degeneracy. It requires the tool bioawk.
It works by looking at the first bases in the read and only prints out the headers that do not have a match with the degenerate primer. Then, comm makes sure that the read was retained in both the R1 and R2 filtering. This allows reads to be filtered whether or not the primers start at the beginning of the R1 or R2 read.
seqtk subseq is where you recover the reads that don't start with the primer.
Be warned though, any off-by-one or sequencing errors will allow some of the sequences with primers at the beginning. Also, this script requires that the headers for the R1 and R2 reads are the same and aren't appended with anything like "\1", "\2", et cetera.
#!/bin/bash
get_primer_free_starts () {
    bioawk -cfastx 'BEGIN{seqs[0] = "CATGGGAGGGCGATCGCAA"; seqs[1] = "CCTGGGAGGGCGATCGCAA"; \
        seqs[2] = "CGTGGGAGGGCGATCGCAA"; seqs[3] = "CTTGGGAGGGCGATCGCAA"} \
        {thisseq = substr($seq, 1, length(seqs[0])); printme=1\
        for (i=0; i<=3; i++) { \
            if (thisseq != seqs[i]){ \
                printme=0 \
            }
        }; if (printme==1){ \
            print($seq)
        }}' "${INPUT} | sort > "${OUTPUT}
}

INPUT="R1.fastq.gz"
OUTPUT="R1_OK.txt"
get_primer_free_starts

INPUT="R2.fastq.gz"
OUTPUT="R2_OK.txt"
get_primer_free_starts

comm -12 R1_OK.txt R2_OK.txt | sort | uniq > passing_reads.txt

seqtk subseq R1.fastq.gz passing_reads.txt > R1.filtered.fq
seqtk subseq R2.fastq.gz passing_reads.txt > R2.filtered.fq


text:
Tool to show DNA sequence and allowing upload of own graph data
Background
We want to be able to load (or request) data for a genome including the sequence and gene annotation (bacteria). Then, we want to load our own annotation which should be displayed as a line plot:
position score
1        5
2        6
etc.

Where each position corresponds to the position in the genome (i.e. FASTA).
Something like this:

Question
Is there a tool that can load a genome sequence, gene annotation, and a position ~ score file, which can be displayed as a line plot?

Welcome to the site. In order to suggest a good alternative could you comment what do you need this tool for? Why do you need to display the gene annotation? What is the problem with IGV ? (I found the sentence a bit confusing, perhaps if you rephrase it I'll understand it)

I thought I was clear regarding what I need it for: uploading genome + annotation + own data (position + score). Why gene annotation? Because we want to look whether the the line plot we upload differs between genes and intergenic regions. Just similar as a GC plot beneath the genes and intergenic regions but then with own supplied scores per position. Indeed, I forgot a word in my sentence about IGV @Llopis

I don't find this question useful. It only helps people who find IGV ugly. Also you don't seem to have researched for other options like IGB. Have a nice day!

A local tool is fine, with upload I meant uploading it to a tool (perhaps not the best word choice). The screenshot is from a tool which only visualizes a Listeria genome (no other bacteria), but I liked the interface. Currently I'm trying to use IGB @terdon

You can use UCSC genome browser, which supports custom genome sequence and annotation. Specifically, one needs to create a Track Hub to hold genomes and annotations.
Here is an example:

http://genome.ucsc.edu/cgi-bin/hgTracks?hubUrl=http://zlab-trackhub.umassmed.edu/yfu/insects/hub.txt&db=hub_129615_cl_v0.8.2

And here is the corresponding hub.txt file holding the metadata, and URLs to actual fasta, bigWig and bigBed files.


text:
How to input data for DESeq2 from individual HTSeq count?
I am comparing the gene expression of 2 bacteria under 1 condition. I have now the count tables for 3 tech. replicates for each bacteria.
Bacteria1_1.count
Bacteria1_2.count
Bacteria1_3.count

...same for the other bacteria.
These files look like this:
gene1 10000
gene2 500
gene3 0
gene4 5000

I want to use DESeq2 for differential gene expression analysis. But I cannot figure out how to properly execute the DESeqDataSetFromHTSeqCount() command with this type of data.
Is there another intermediate step to add?

What errors are you seeing? How are you trying to do this? Please [edit] your question to show us the commands you tried.

Probably you should join all the count tables into a single one (If they share the same rows)

cross posted here

In case still needed, here is a piece of code that I have used for the same purpose (includes a lot of sanity checks). This is highly likely an overkill as each count file would probably have the same rows (as a result of being generated from the same workflow) and thus can easily be merged column-wise.
count_file_dir <- "../input/count_files/"
count_file_names <- list.files(count_file_dir)
count_files <- paste0(count_file_dir, count_file_names)

# read and merge count files, while removing last 5 lines that correspond to summary stats
count_table_list <- mclapply(count_files, function(x) {fread(x)[1:(.N-5)]})
names(count_table_list) <- sub(".count", "", count_file_names)

# check the dimensions of count files
t(sapply(count_table_list, dim))

# check if the count files have the same gene ordering
# compare gene names of the first file with all
gene_names <- sapply(count_table_list, function(x) x$V1)
as.data.frame(apply(gene_names, 2, function(x) all(x == gene_names[,1])))

# add file name as a unique column name for gene counts
for(i in names(count_table_list)) {
  names(count_table_list[[i]])[2] <- i
}

# merge read counts and convert to a matrix
count_matrix <- Reduce(merge, count_table_list)
count_matrix_rownames <- count_matrix$V1
count_matrix$V1 <- NULL
count_matrix <- as.matrix(count_matrix)
rownames(count_matrix) <- count_matrix_rownames
colnames(count_matrix) <- sub(".count", "", colnames(count_matrix))

dim(count_matrix)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = ...,
                              design = ...)


text:
Machine learning using protein-sequences
I'm participating in a bioinformatics machine-learning seminar at my university. The main task is predicting binary classification of protein-protein interactions using sequence data as input.
One of the subtasks is familiarization with the dataset and presenting the dataset. Now I'm wondering which additional information I can get out of the sequence data.
I want to start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact. The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”. Species are completely mixed.
Is there something interesting that I can additionally include in my analysis? I'm used to exploratory data analysis though this is different due to the nature of the protein sequences.

Hi there, it would be perhaps good to describe the problem bit more. Do I get it right that you have a dataset of amino acid chains and now you want to predict witch do proteins interact? Do you have a gene position of proteins in the genome? Expression data? Do you have species? What is your training dataset?
These might be completely irrelevant questions since I am not a proteomic person. However, I don't really understand the question as it is right now.

Hi Kamil, thanks for the follow up question!
Yup, you get the problem right! I wanna start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact.
The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”.
Species are completely mixed.

Great, we like juicy details. Could you perhaps [edit] your question and add it there?

Any specific kind of machine learning algorithm you intend or have to use?

Hi Florian, for our assignment we use ANN as well as SVM

There are two types of features which can be extracted from protein sequences.  It is not possible to known which, if any, of these features would be useful in classification tasks.  It may not be possible to build a classifier at all, or it may be very straight forward.  To know this, a feature selection technique must be used, such as e.g. forward feature selection.
Sequence based features
Sequence based features are derived from sequence-only information.
The most common used features are the PseAAC (pseudo amino acid composition) style features, including:

Average Position of each amino acid - 20 features
Average Distance between each amino acid - 20 features
Motifs - 20 features, 20*20 features, and/or 20*20*20 features (depending on motif length)
OAAC (Overall Amino Acid Composition) - 20 features
SAAC (Split Amino Acid Composition) - 60+ features
AAIndex Composition - 1 to 567+ features
PCC (Physicochemical Composition) - ~13 features
PSSM encoded as 1, 20, 200, 400 or 600 features (i.e. Psi-Blast output using NR or swissprot)
Average predicted residue burial/exposure level
Average predicted sequence entropy

Structure based features (predicted)
Using the sequence, you can either predict the 3d structure using homology, ab initio, or directly use a homologous structure.  You can then extract predicted 3d features, such as:

3d size (i.e. max of x,y,z)
Calculated ASA/SASA
3d contacts (type, frequency, etc)
any other local/global environmental information
cath/scop classification
dssp/stride composition


text:
Generating basic dna sequences in R
Looking for some guidance.
In an introductory bioinformatics class I am tasked with doing the following in R:

Create the sequence ACTGACTG
Generate a random string of nucleotides with equal probabilities.

Are there certain functions that I should be looking at? No library or package or set of functions was provided so I am completely lost. I am assuming there must be functions in R that do these things that I need to manipulate, but I cannot seem to find any resources explaining that. I looked at seq() (thinking it might generate DNA sequences), but it only seems to generate sequences of integers. I also looked at randDNA(), but it is random, how could I make it generate the exact string in question 1).
Any direction will be greatly appreciated.
Thanks.

They are two separate problems.  If you have something that generates sequences of random integers can you think of a way to map those to the bases you want to report.  The first question is ambiguous.  How do your report the sequence you created?

Without giving away too much for your homework check out the sample() function in R. Think about how you can pass DNA bases to this function to generate random sequences.

Thanks this is really helpful. Can you please help me understand the prob parameter, as I understand the probability of each element is relative to the entire list. Is it possible to specify probabilities for each element in certain positions? say I have sample(1:2, size=4, replace=TRUE, prob=c(0.8, 0.2)) This says that 1 has probability 0.8 is the list overall, but I want to say 1 has probability 0.8 in position 4 (and some other probability in other positions)

This can't be done using just the prob parameter, because it only takes a single vector as input. You'll have to find some loopy, or otherwize crazy way for trying to apply different probabilities at different positions.


text:
Store contigs and genome substrings based on mismatches
I have a string, a list of substrings, and a list of positions where substrings are found in string.
string = 'ATGGTACACGCTACGAGCTAG'
substrings = ['TACA', 'CTAGAAAAG']
positions = ['11, 4', '']

Substring 'TACA' appears at positions 11 and 4 in string (i.e. '11, 4').
Substring 'CTAGAAAAG' does not appear in string (i.e. '').
I need to remove 'CTAGAAAAG' from substrings and I need to remove '' from positions, while creating a new list with 'CTAGAAAAG'. I need your help here to create code which will give the following results:
new_substrings = ['TACA']
new_positions = ['11, 4']
new_list = ['CTAGAAAAG']
string = 'ATGGTACACGCTACGAGCTAG' (obtained from above)

Once this is complete, the given code:
new_positions = [map(int, x.split(',')) for x in new_positions]
seq_pos = dict(zip(new_substrings,new_positions))
seq_ref = {k: [string[start:start+length] for start in starts] for k, starts, length in zip(new_substrings,new_positions,[len(x) for x in new_substrings])}

Gives:
seq_pos = {'TACA': [11, 4]}
seq_ref = {'TACA': ['TACG', 'TACA']}

Then, from seq_ref, I need to create a list of integers which represents the number of mismatches between the key: 'TACA' and both of its values: 'TACG' and 'TACA'.
mismatches = [[1, 0]]

The code to find the # of mismatches between two strings is provided:
def HammingDistance(p, q):
    dist = 0
    for i in range(len(p)):
        if p[i] != q[i]:
            dist += 1
    return dist

Finally, I need to save my data in a CSV file (i.e. Excel or Numbers).
new_positions = ['11, 4'] (obtained above)
new_substrings = ['TACA'] (obtained above)
actual_substrings = [['TACG', 'TACA']] (values of seq_ref)
mismatches = [[1, 0]] (obtained above)

Info needs to be sorted based on # of mismatches.

I understand this is a multi-step problem but any help would be appreciated.

You posted a lot of given code, but what have you tried so far? To be honest this reads a lot like a homework assignment and while most users here are happy to help this is not a site to get your work done for you.

Where did the blocks of "given" code come from?  I agree this seems like homework and many of the questions you are trying to answer have similar concepts in the "given" code.  Have you tried using or understand some of the functions in the "given" code to see how they apply to the missing steps you have?

This is not part of a homework assignment.  I've never taken a coding or bioinformatics class before and I just need some assistance.  This is a unique problem that I have created from scratch.  The given code is my own work and I uploaded it so that it may provide a better understanding of the problem at hand.  I need help removing a blank string element from a list and I need help iterating the HammingDistance() function for elements in my list.  Also, I am unfamiliar with csv files, so I am just reaching out to see if anyone can help.

Ok I posted a solution.  It is confusing because most of what you are stuck on short of the CSV are concepts that you are using in the "given" code, including zip, dictionaries, list comprehension / iteration.

For the first question try something like:
new_substrings = []
new_positions = []
new_list = []
for pos, sub in zip(positions, substrings):
   if not pos:
       new_list.append(sub)
   else:
       new_substrings.append(sub)
       new_positions.append(pos)

this gives:
>>> new_substrings
['TACA']
>>> new_positions
['11, 4']
>>> new_list
['CTAGAAAAG']

For the hamming distance you just need to iterate over the values in seq_ref and append the results of the hamming distance function in mismatches.  Based on the given code and the example above you should be able to create that block.  Below is an untested quick example to get you on track.
mismatches = []
for seq in seq_ref.keys():
   for value in seq_ref[seq]:
       mismatches.append(HammingDistance(seq, value))

As far as the CSV and sorting is concerned I suggest that you look into pandas or numpy.  Both of those packages provided functions to help with this.

Just a note to say that your hamming distance calculation is pretty slow. When we've had to do this calculation before on many millions of sequences we've ended up cythonising the code. But friends who are better coders have even better algorithms for calculating these things.


text:
Error while performing "lift-over" of genomic coordinates between two species using Ensembl Compara API
I'm trying to convert coordinates between the reference mouse genome and CAST/EiJ mice strain. Since Ensebml already has alignment available in the browser, I decided to use Ensembl perl API to lift over the alignment for the regions I'm interested in. The problem is that it doesn't seem to be working. I'm getting following error message:
-------------------- EXCEPTION --------------------
MSG: The database do not contain any CACTUS_HAL data for Bio::EnsEMBL::Compara::GenomeDB=HASH(0x557db0d661e8)!
STACK toplevel coord.pl:42
Date (localtime)    = Mon Oct 15 07:42:17 2018
Ensembl API version = 94
---------------------------------------------------

which I assume means that I'm either searching for the wrong alignment record or the alignment itself is not available. Following is the perl script which I modified from the Ensembl Compara API tutorial:
use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Exception qw(throw);
use Bio::AlignIO;
use Getopt::Long;

my $species = "murinae";
my $targetSpecies = "mus_musculus_casteij";
my $alignment_type = "CACTUS_HAL";
my $output_file = undef;
my $coord_system = "chromosome";
my $seq_region = "14";
my $seq_region_start = 5000000;
my $seq_region_end = 5010000;
my $output_format = "clustalw";
my $genome_dbs;

use Bio::EnsEMBL::Registry;
Bio::EnsEMBL::Registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org',
    -user => 'anonymous',
    -port => 5306);

my $genome_db_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomeDB');

throw("Cannot connect to Compara") if (!$genome_db_adaptor);

my $this_genome_db = $genome_db_adaptor->fetch_by_name_assembly($species);
my $target_genome_db = $genome_db_adaptor->fetch_by_name_assembly($targetSpecies);
push(@$genome_dbs, $this_genome_db);
push(@$genome_dbs, $target_genome_db);

my $method_link_species_set_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'MethodLinkSpeciesSet');

my $method_link_species_set =
$method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
      $alignment_type,
  $genome_dbs);

throw("The database do not contain any $alignment_type data for $target_genome_db!")
    if (!$method_link_species_set);

my $slice_adaptor = Bio::EnsEMBL::Registry->get_adaptor($species, 'core', 'Slice');

throw("Registry configuration file has no data for connecting to <$species>")
if (!$slice_adaptor);

my $query_slice = $slice_adaptor->fetch_by_region('toplevel', $seq_region, $seq_region_start, $seq_region_end);

throw("No Slice can be created with coordinates $seq_region:$seq_region_start-".
    "$seq_region_end") if (!$query_slice);

# Fetching all the GenomicAlignBlock corresponding to this Slice:
my $genomic_align_block_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomicAlignBlock');

my $genomic_align_blocks =
$genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice(
      $method_link_species_set,
  $query_slice);

my $all_aligns;

# Get a Bio::SimpleAlign object from every GenomicAlignBlock
foreach my $this_genomic_align_block (@$genomic_align_blocks) {
    my $simple_align = $this_genomic_align_block->get_SimpleAlign;
    push(@$all_aligns, $simple_align);
}

# print all the genomic alignments using a Bio::AlignIO object
my $alignIO = Bio::AlignIO->newFh(
-interleaved => 0,
-fh => \*STDOUT,
-format => $output_format,
    -idlength => 10
);

foreach my $this_align (@$all_aligns) {
    print $alignIO $this_align;
}

exit;

Any pointers on what I'm doing wrong ?

Perl is a language I haven't heard in a long time :D It seems like the GenomeDB is passing an argument that might be wrong. Are you using the latest bioperl/perl?

While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks.
The quick fix - in python looks like this:
#This is the main function which uses Ensembl REST API to access the Enesembl Compara database
#And retrive the alignment blocks which correspond to the MM10 standard.
#The function accepts chromosome, start and array and returns the alignment blocks
#reported by the API. Since the blocks can be on different chromosomes/and/or/contigs, it
#returns the whole output in form of dataframe - containing chromosome, start and stop of the
#alignment blocks corresponding to input and all the other alignment blocks corresponding the
#16 mice species, if available.
#!!!!!!!----------Caution start----------!!!!!!!
#The API is a bit buggy, so it is advised not to request for specific species.
#!!!!!!!-----------Caution end-----------!!!!!!!
def liftover(chrom, start, end):
    mmStart, mmEnd, mmchr = [], [], []
    caStart, caEnd, cachr = [], [], []
    #ryStart, ryEnd, rychr = [], [], []
    server = "https://rest.ensembl.org"
    ext = "/alignment/region/mus_musculus/" + str(chrom) + ":" + str(start) + "-" + str(end) + ":1?species_set_group=murinae;method=CACTUS_HAL;content-type=application/json"
    headers={ "Content-Type" : "application/json"}
    r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
    if(r.ok):
        json_string = json.dumps(r.json())
        jData = json.loads(r.content)
        results = r.json()
        print("The response contains {0} properties".format(len(jData)))
        print("\n")
    else:
        r.raise_for_status()

    for element in jData:
        element.pop('tree', None)

    for element in jData:
        for key, value in element.items():
            for z in value:
                if (z["species"]=="mus_musculus"):
                    mmStart.append(z["start"])
                    mmEnd.append(z["end"])
                    mmchr.append(z["seq_region"])
                if (z["species"]=="mus_musculus_casteij"):
                    caStart.append(z["start"])
                    caEnd.append(z["end"])
                    cachr.append(z["seq_region"])
    if (len(mmStart) != len(caStart)):
        print("Length mismatch between species blocks - proceed with caution!")
    return [min(mmStart), max(mmEnd), min(caStart), max(caEnd), len(mmStart)]

Again, it depends on how alignment blocks are reported - and requires further digging to resolve which mm10 alignment block resembles CAST alignment block.


text:
What is a site in a phylogenetic tree?
When I read about a term related to phylogenetics, for example, $K_a$/$K_s$ ratio, I constantly see the word site, whose meaning isn't explained and it puzzles me.
For example substitutions per synonymous site, substitutions per non-synonymous site, substitutions at a site etc.
Could you define a phylogenetic site?

So there were 0.6 substitutions per the blue-violet site, what does it mean? What is the length of the site? Is the length of sites constant? Is it the length of a sequence read? Are introns the part of the site?

A site is simply an individual discrete position, normally a single nucleotide (or amino acid). For example, consider a toy alignment:
sequence1 ATGC
sequence2 ATCC
#           ^ arrow points to the third site

Here there is 1 mutation in 4 sites, so 0.25 substitutions per nucleotide site.
I'm not sure where the term originated, but probably one of the earlier definitions is from renowned population geneticist Motoo Kimura, who simply stated in 1968:

In this paper, “site” refers to a single nucleotide pair

A pair in this context, means either a homozygous or heterozygous site.
Synonymous sites are positions where nucleotide changes lead to no amino acid change, these are typically at the 3rd position of the codon where the nucleotide is often degenerate, for example, TAT, TAC TAG, and TAA all encode the Serine amino acid. By contrast non-synonymous sites are where a single nucleotide change does lead to an amino acid change, and therefore possibly a more important functional change to the resulting protein.
Thus, the ratio of Ka/Ks (or Dn/Ds) is thus often used as a crude proxy of selection, with Ka/Ks < 0, Ka/Ks = 0 and Ka/Ks > 0 implying purifying (negative) selection, neutral evolution and positive selection respectively. Note that this assumes the neutrality of synonymous sites, which is increasingly questioned due codon usage bias and regulatory elements overlapping coding regions (and also ignores mutational biases and biased gene conversion, although these can be accounted for in an evolutionary model, and also assumes no sequencing errors and perfect alignments, which can impact a lot on the results too).


text:
vcftools: indel size histogram command returns empty file
I would like to get some summary statistics on a vcf file from one individual, which has over a million variant calls. I've tried to make a histogram of indel sizes with this command,
vcftools --vcf sample.vcf --out sample --hist-indel-len

but it returns an empty file. I see thousands of deletions in the vcf, although strangely no insertions. A typical line in the vcf file I have looks like
NC_025968.3 168363 191 N <DEL> . . SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0 GT:SU:PE:SR ./.:19:19:0

Update: I aligned reads to my reference with BWA MEM, then used samblaster to extract discordant pairs and split-reads. I sorted and indexed the three bam files, then used lumpyexpress to generate the vcf.
samblaster -i sample.sam --excludeDups --addMateTags --maxSplitCount 2 --minNonOverlap 20 | samtools view -S -b - > sample.bam

samtools view -h sample.bam | samblaster -a -e -d sample.disc.sam -s sample.split.sam -o /dev/null

lumpyexpress -B ../sample.sorted.bam -S ../sample.split.sorted.bam -D ../sample.disc.sorted.bam -o sample.vcf

I would like to see the size distribution of deletions, and if possible, figure out why there are no insertions.

Thanks for your question, Joanne, and welcome to Bioinformatics SE. Can you post some lines from your VCF file so that we can get a better idea of what is going on? It is hard to tell without looking at the data directly.

Thanks for the welcome! Here is a typical line: NC_025968.3     168363  191     N       <DEL>   .       .       SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0        GT:SU:PE:SR     ./.:19:19:0

Hello @Joanne just a note. You can [edit] your question to add details (like the line I have already added for you) or any other progress or idea you will have regarding the question. Usually more details you give, better chances that there will be somebody who will know the answer.

Just curious, how are you calling this VCF file? If you are using GATK I think there are some options to output only SNPS, or SNPs + indels, or just indels. Also, can you clarify exactly what you mean by summary statistics? Right now it is hard to tell what your question is exactly: IE do you want a size distribution of the deletions that are in there, or do you also want to figure out why insertions are missing from your VCF? If you would like the latter, we'll need info on the commands you used to generate the VCF file.

Thanks for the forum protocol tips - I have edited my question to include more details.

It looks like lumpy looks for breakpoints in your reference. Presumably places where the is a haplotype switch detectable with the data on-hand?
My only thought is that since lumpy is being used to produce the VCF file that it is inherently not looking for insertions. I would suggest running the same BAM file through the GATK best practices pipeline (easier said than done) and vcftools pipeline. If those produce insertions but lumpy does not, it must be due to how lumpy produces the calls. Good luck!


text:
How to assign the best gap penalty and gap extension penalty using BLOSUM65
For an assignment I must do a pairwise optimal local alignment using BLOSUM65 and five protein sequences. The algorithm I want to use is the Smith-Waterman.
Context protein sequencing using Blastp: all the sequences are the hemoglobulin subunit a or b;

two times two pairs are from the same organism (but differ in subunit),
one is not;
all the sequences differ in sequence length (around ~15 amino acids. Gaps!)

The default for the BLOSUM65 penalties are: 10 for a gap, 1 for an extension.
How do I optimise the gap penalty and gap extension penalty? Should I stick with the default? If so, why?

I'm voting to close this question because it doesn't seem to have enough detail. Smith-Waterman is a well-known algorithm, the question is explicitly a homework question, but it's not clear from the question what form the answer should be in (e.g. what programming language). There's confusion because the question introduces BLASTp, which is a separate algorithm from Smith-Waterman.


text:
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).

Welcome to the group.  Can you expand on your question to include what tools you are using, what you have tried.  How did you generate the curve and do you have any sample data to help answer your question?

Hi Bioathlete. I'm totally new in Bioinformatics field. I have not tried much yet due to lack of knowledge, but I edited my question. Thank you.

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.


text:
Script to allow gene set enrichment analysis of 10x genomics data in R
I have 10x single cell RNA seq data. Which R package is best suited for analysis of the 10x data matrix. What is the script to prepare the data for downstream GSEA analysis.
I have already processed samples for single cell data (for multiple other purposes). my intention is to focus on subpopulations i.e. cell type x at different time points i.e. time 0, 1, 2 etc and compare the enriched gene pathways in this one cell type at different time points. there are three files that are used (at least in the R Seurat package) with the following file extensions: .mtx, .tsv (barcodes.tsv, genes.tsv and matrix.mtx). the R Seurat package merges these files to create a "Seurat object" which helps mainly with identifying individual diffl gene expression in each barcode cell
I found a package for R: goseq that at least has been used for bulkRNA seq but I was hoping someone has used this package specifically for single cell RNA seq data, and what script has been used to merge the 3 files to create an object for analysis. pardon my terminology, but I am a novice to R and have experience primarily through package use.

Welcome to the group.  Can you specify what format you data is in and what format you need it in for the downstream GSEA.  What tools have you tried or looked into?

Hi Jay, its not clear from your question what the biological question you are trying to address with a GSEA analysis. Do you want do do enrichment of the genes that are expresesd on average in the cells of this experiement? If this is the case, how would this differ from an analysis of bulk sequencing? Do you wish to know the pathways that differ between different scRNAseq samples? Do you wish to know the pathways that differ between clusters of cells within one sample? The biological question at hand will strongly influence the downstream analysis approach.

In addition to Ian question, you mention goseq. Are you aware of the differences between pathways and gene ontologies terms? Do you have in mind any pathway database or you don't care (there are several with different focus). Are we talking about human pathways or other organisms? For other organisms there are fewer databases of pathways or less consensus. Keep updating the question. I think it can become a great question for future readers

Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.
Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable.  The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.

There is no purpose-built R package to perform gene set enrichment analysis on single-cell data but there does not need to be. You should be able to tools developed for bulk-RNA-Seq or microarray data, although you may not get as significant results from a sparse scRNA-Seq matrix as single-cell technologies have poor sensitivity and miss genes. What you need for a gene set enrichment analysis is:

a database of pathways (or gene sets) to test enrichment for. Packages such as GO.db or reactome.db from Bioconductor support these in R.
a gene expression matrix which you can get from the output of cellranger (as matrix.mtx). This should be normalised and corrected for batch effects. The Seurat R package is recommended to do this. Calling object@data will return the normalised seurat data matrix to pass it to other tools.
an annotation to identify the genes (depending on the reference data you used for the cellranger "transcriptome").
to define groups for differential expression analysis. You can import the clusters generated by cellranger in outs/analysis or export clusters of your choosing from the Loupe Browser as CSV files that can be imported into R. These should match the barcodes.tsv of columns of your data matrix to define the cluster for each cell barcode.
to run a differential expression analysis to compare clusters. You can use clusters to define cell groups and perform analysis with tools you would use for bulk-RNA-Seq such as limma, edgeR, or deseq. Don't forget to adjust for multiple tests. You also need to correct for batch effects between replicates and log-transform you data to ensure the results reflect biological differences. This will give you a list of differentially expressed genes.
to calculate gene set enrichment of differentially expressed genes in a gene set list such as those supported by pathway databases. You can do this like any other list of significant genes with tools used for bulk-RNA-Seq data, either with dedicated webtools or R packages such as safe or GSEA. You can also run the hypergeometric test in R yourself.

The cellranger output gives you all the data you need to use standard tools for which there are already guides on how to do differential analysis. These analyses do not differ from bulk-RNA-Seq analysis once the data has been normalised and the  clusters have been defined.
[edit:corrected typo in package name]

I realized that your question is specifically gene set enrichment analysis. However, I personally have never tried any R package doing this. Although topGO may seem to be the one that may work for you? What I have been doing is that for gene enrichment analysis using Gene Ontology, simply paste and run, or DAVID when it was still updating. Sorry, I am not an expert in incorporating biological significance study in the computational result.
It seems that you are also not quite sure what tools you should use for downstream analysis. I think the following publication may be helpful.
https://www.ncbi.nlm.nih.gov/pubmed/30228881
The same primary author also published a paper titled, Cluster Headache: Comparing Clustering Tools for 10X Single Cell Sequencing Data (For some reason I could not find the link to it, will edit the answer if I find it. )
Hope that helps.

DAVID is widely regarded as outdated. The databases that it uses are not kept up to date. This is not recommended.

@TomKelly That is absolutely right. That's why I said it was good when it was updated. Not recommended until it starts to update again.


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

Welcome.  What tools have you been using?

hubalign ......

Do you need it in python? I know one very nice package in R that does that.

@Llopis it would be nice in python but if any R package you know, kindly mention that

You can calculate gene ontology similarity with the GOSemSim R package (paper)(disclaimer I'm a contributor to this package). You have several similarity scores implemented, some of them are across the three subontologies.
In python you have the GOATOOLS (paper). I think there was another Python package but now I don't remember its name.

thanks @Llopis. i looked into goatools already but i don't think that they are have methods to calculate GO CC, MF or BP

Please, mention in the questions what you have done, because now I feel that I wasted my time trying to answer you when your problem is something else. Did you look at the notebook ?

sorry if you feel that way. but did you saw any GO calculated in goatools ?

I linked a notebook explaining how they do it..


text:
Does read sequence identifier need to be unique for proper alignment?
I am preprocessing a set of scTHSSeq reads in which each read sequence identifier has been replaced with the cell barcode, e.g.:
@CGAGGCTG_TATCCTCT_AGTCCA
...
...
...
@CGAGGCTG_TATCCTCT_GGTACG
...

The read sequence identifier has been replaced with the cell barcode to be able to demultiplex alignments after bwa into a different bamfile for each cell.
I ran the alignment using bwa mem and I received no errors. However I would like to know if there are caveats or reasons why this approach may lead to misleading results.
I will then proceed to remove PCR duplicates with samtools rmdup executed in each cell's bamfile and then I'd run peak calling using spp.

What tool did you use to replace the read names? Can you modify that to split up the reads into separate files or encode the cell barcode in the read description instead of the read name?
Typically is it assumed by many tools that the FASTQ read names are unique. you could have issues with BWA and split alignments but I don't know if that pertains to your data.  I think the best you can do is the sanity check the results at each step of the process.
In the fastq format (https://en.wikipedia.org/wiki/FASTQ_format) on the first line denoted by the '@' there can be two values separated by a space.  The first is the name and the second is the description.
So an example would be:
@Original_read_name CGAGGCTG_TATCCTCT_AGTCCA
ATCGATCATGATGACA
+
@&%*@(%@%*($!##!

In this iteration I did it using awk, do you mean incorporating it in the first line after the sequence identifier? I can certainly modify the script.

I mean the one defined as entry.comment here?

I edited my answer to address your question

I noticed now that BWA does not preserve read description, I should probably append cell barcode to read name or this may introduce problems?

The classic way to do this would be to append the barcodes to the read name such that the read name looks like: @Original_read_name_CGAGGCTG_TATCCTCT_AGTCCA, If the barcodes are extacted from the read sequences themselves, this processing can be achevied using umi_tools extract


text:
Converting Gene Symbol to Ensembl ID in R
I'm trying to convert ~20,000 different human gene symbols to ensembl IDs. I've been trying to use biomaRt to do this, but continue getting the following error
getBM( attributes=c("ensembl_gene_id") , filters= "mgi_symbol"  ,mart=ensembl)
Error in martCheck(mart) :
  No dataset selected, please select a dataset first.  You can see the available datasets by using the listDatasets function see ?listDatasets for more information.  Then you should create the Mart object by using the useMart function.  See ?useMart for more information

Is there a different format I should be using for this? I tried using BioConductor's description of how to use getBM but am a R novice and confused as to how to set this up to convert gene symbols into ensemble IDs.

You say you have human gene symbols in your text, but in your code you say you are filtering by MGI symbols. Which is it?

How can I do if I have gene accesion number as NM_? thank you very much

You can use the standard annotation package for humans (you can also use biomaRt, but it can be more confusing, see below):
library("org.Hs.eg.db") # remember to install it if you don't have it already
symbols <- mapIds(org.Hs.eg.db, keys = ensemblsIDS, keytype = "ENSEMBL", column="SYMBOL")

The ensemblsIDS should be a character vector: c("ENSG001", "ENSG002", ..).
However, if your primary key is an ENSEMBL id it is preferred to use a ENSEMBL package centred around the ENSEMBL id :
library("EnsDb.Hsapiens.v86")
mapIds <- mapIds(EnsDb.Hsapiens.v86, keys = ensemblsIDS, keytype = "GENEID", columns = c("SYMBOL")

The solution with biomaRt:
library("biomaRt")
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
getBM(attributes='hgnc_symbol',
      filters = 'ensembl_gene_id',
      values = ensemblsIDS,
      mart = ensembl)

Thanks much for the detailed response!! If I currently have ensembleIDs as a .csv file in columns, what is the best way to convert them to a character vector? Also, just to confirm, ensemblsIDs (defined as the value) is the input data set, correct?

Well, read around the forums. Normally one would read the .csv file and then extract the values with something like as.character(csvData[, "ensemblColum"]). It should be correct, unless I understood something different.

I tried several R packages (mygene, org.Hs.eg.db, biomaRt, EnsDb.Hsapiens.v79) to convert Ensembl.gene to gene.symbol, and found that the EnsDb.Hsapiens.v79 package / gene database provides the best conversion quality (in terms of being able to convert most of Ensembl.gene to gene.symbol).
# Install the package if you have not installed by running this command:
BiocManager::install("EnsDb.Hsapiens.v79")

library(EnsDb.Hsapiens.v79)

# 1. Convert from ensembl.gene to gene.symbol
ensembl.genes <- c("ENSG00000150676", "ENSG00000099308", "ENSG00000142676", "ENSG00000180776", "ENSG00000108848", "ENSG00000277370", "ENSG00000103811", "ENSG00000101473")

geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v79, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))

# 2. Convert from gene.symbol to ensembl.gene
geneSymbols <-  c('DDX26B','CCDC83',  'MAST3', 'RPL11', 'ZDHHC20',  'LUC7L3',  'SNORD49A',  'CTSH', 'ACOT8')

geneIDs2 <- ensembldb::select(EnsDb.Hsapiens.v79, keys= geneSymbols, keytype = "SYMBOL", columns = c("SYMBOL","GENEID"))

Other available R packages / gene database for the conversion can refer to this GitHub page.


text:
What are "split reads" and "intron clusters?"
I'm working through the example data set for LeafCutter and the documentation mentions "split reads":

This will cluster together the introns fond [sic] in the junc files listed
  in test_juncfiles.txt, requiring 50 split reads supporting each
  cluster and allowing introns of up to 500kb.

Searching "split reads" on Google yielded few results, and the results I did find were hard for me understand,

Chimeric reads occur when one sequencing read aligns to two distinct
  portions of the genome with little or no overlap. Chimeric reads are
  indicative of structural variation. Chimeric reads are also called
  split reads.

One sequencing read aligns to two distinct portions of the genome, but then with little or no overlap? Isn't that a self-contradiction?
I can more or less deduce what "intron clusters" are but I figured I might as well try to solicit a more formal definition. I know what introns are and I can figure what it means to cluster them together. I guess I just want to know what split reads and intron clusters mean in the context of working out the example data for LeafCutter.

split reads - These are read that have two or more alignments to the reference from unique region of the read.  In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' end of exon2 and bases 76-150 aligning to the 5' end of exon3.  This would be a split read because it have two alignments (exon2 and exon3) and those alignment are from unique (non overlapping) parts of the read, 1-75 and 76-150.
Chimeric reads are an example of split reads.  Many people use "chimeric reads" to describe a type of noise that occur during PCR or sequencing where the bases from two different molecules are reported as one read and the alignments of that read become split when it aligns back to the location of the two original molecules.
"intron clusters" - In this case I think that the LeafCutter software is describing stacks of split reads to cluster to support the detection of an intron.

"Chimeric reads are and example of chimeric reads" - there are, um, at least two problems with that sentence!

hey, so can we say that a split read can be a candidate of a gene/transcript fusion?

Yes split reads could be.  Many fusion detection tools use split read information


text:
Searching for start and stop codons for protein sequencing of contigs
I need to convert contigs into their respective protein sequences given a reference genome (i.e. I need to take a substring, whose position is already known on the string, and I need to locate the nearest start and stop codons).
This is tricky because sometimes the first position of triplet codon within the contig will not be a multiple of three and could be located in the intergenic region (i.e. between genes). The goal is to separate both coding and non-coding DNA.
This is my code so far:
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna, generic_protein

start_codons = ['ATG']
stop_codons = ['TAG', 'TAA', 'TGA']

string = 'GG*TAG*CCAATT*ATG*AACGAA*TAG*GAC' #remove '*', just for visual
contigs = ['CCAA', 'TGAAC', 'GAA', 'GGAC']
positions = [5, 12, 17, 22] #position indices for each contig on string

extended_contigs = []
extended_position_contigs = []
intergenic_contigs = []
intergenic_position_contigs = []
for i in contigs:
    extended_contigs.append(#some code)
    extended_positions_contigs.append(#some code)
    intergenic_contigs.append(#some code)
    intergenic_positions_contigs.append(#some code)

I should get extended_contigs = ['ATGAAC', 'ATGAACGAA'] and extended_positions_contigs = [12, 17].  These are the contigs that are located within genes. In order to code them into peptides, I need to reach back in the string until I find the start codon and extend the initial contig (e.g. TGAAC -> ATGAAC and GAA -> ATGAACGAA)
I should also get intergenic_contigs = ['CCAA', 'GGAC'] and intergenic_positions_contigs = [5, 22].  When the missing code is run, the computer searches to the left of the string and finds a stop codon (e.g. TAG) before the start codon. Thus, the contig is located in between two genes and nothing needs to be added. These intergenic contigs are just stored in a new list.
My code continues:
prot_contigs = []
for i in extended_contigs:
    my_dna = Seq(i, generic_dna)
    my_prot = my_dna.translate()
    prot_contigs.append(str(my_prot))

Here, no new code needs to be added. After the above is run, prot_contigs = ['MN', 'MNE'].
The last step of code (which I need help with) will transform prot_contigs into new_prot_contigs = ['MN', 'E'].
How? If for any contig (e.g. 'TGAAC'), the beginning or end is part of another codon (not perfect multiple of 3), the extra codon on either end will be maintained (e.g. 'MN' stays 'MN').  Otherwise, if the contig (e.g. 'GAA') perfectly matches a codon, anything that was added onto it will be removed (e.g. 'MNE' becomes 'E').
I would try to solve the 2 portions of code myself, but I am unsure on how to take a position on a string (i.e. beginning points of contig) and look along the string to find the nearest start/stop codons, so I could determine the DNA's function and accurately sequence the protein coding contigs into peptides.
Any help would be appreciated!

Are you trying to predict open reading frames and you don’t want to use someone else’s gene prediction tool?

@Chris_Rands Unfortunately, I am not familiar with open reading frames or gene prediction tools.  Would such a software be able to solve the above? If so, is there open source code I could use?

@ThomasRibeiro you can explore the code of TransDecoder which is pretty widely used.

If your string is 'GGTAGCCAATTATGAACGAATAGGAC' then I would assume that you could separate this into coding and non-coding regions such as:
GG*TAG*CCAATT*    (non-coding)
ATG*AACGAA*TAG    (coding region)
*GAC              (non-coding)

It starts with a start codon ATG and ends with an end codon TAG.
Your first event position 12 does not make sense to me. What does * mean? It seems of that your sequence starts with GG which is not a triplet? Perhaps did you mean GGG?
I'm the author of node based DNA visualisation AminoSee. It has a bug in that it only renders half right: because DNA can be read backwards, and my program only does it forwards. I actually came to BioInformatics today to try and get my head around reading backwards and how that might work.
AminoSee only transcodes 3 base pairs at a time, and discards triplets that aren't, well, 3 base pairs longs. I'm self taught and so I guess unaware if it's possible to have any sequence that is not a multiple of 3. Perhaps it's possible. I assumed not. My program is more of an art piece that a science tool, but I plan to use it one day to perform a "visual diff" between two species (Humans and Chimp). See AminoSee  DNA Visualisation over at aminosee.funk.nz


text:
How to convert featureCounts to FPKM?
I have seen many posts regarding counts to RPKM and TPM. I haven't seen any post for counts to FPKM.
I have RNA-Seq data which is paired-end reads. Extracted the counts using featureCounts for all the samples.
There is a function to convert counts to RPKM: using the gene_length
rpkm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(counts) * 1e6
}

I know that RPKM is mainly used for single-end reads data. Do you think I can use the above function for converting counts to FPKM as my data is paired end? [TCGA data: HTSeq counts are converted to HTSeq FPKM]
Somewhere I have seen that the same function applies for single end data which will be RPKM and also for paired end data which will be FPKM. Am Ir right or wrong?
If not can anyone show some function or code to convert counts to FPKM please.
thanq

Why do you need RPKM/FPKM? Haven't you heard that these values are obsolete?

I guess OP wants to follow the most up to date procedures, but seeing that in TCGA they used FPKM is trying to work like them. I thought that FPKM were ok: see this https://bioinformatics.stackexchange.com/a/4270/48, or this https://bioinformatics.stackexchange.com/a/69/48. @beginner I hope the links are helpful, update the question otherwise

@Llopis thanks for the links. In this https://bioinformatics.stackexchange.com/questions/66/how-to-compute-rpkm-in-r/69#69 I see that FPKM function using effective length...it is gene_length right?

@b.nota Yes, but as TCGA data is FPKM I wanted to use FPKM data.

@beginner See this question and answers: https://bioinformatics.stackexchange.com/q/367/48

I have seen many posts regarding counts to RPKM and TPM.

There’s your answer then: FPKM = RPKM. It’s simply a more accurate name.
Speaking of RPKM for paired-end data is discouraged because the reference to “read” in this context lends itself to ambiguity. But mathematically the quantity is the same: we are counting fragments, not individual reads (of which each fragment has two, for paired-end data).
But as mentioned in the comments, there are good reasons against using FPKMs.

thanks for the answer. I have a question. I have Ensembl gene ids as rows and samples as columns with counts. I have an other column gene_length. I'm trying to convert those counts to TPM. In this post [https://gist.github.com/slowkow/6e34ccb4d1311b8fe62e] they used gene_length to convert counts to TPM. But I have seen one of your post in Rpubs [http://rpubs.com/klmr/rnaseq-norm] where you used effective length. Which one is the right one for converting counts to TPM?

In every case, the correct number is the effective length, i.e. the transcript/gene length minus the average read length plus 1 (this is the “effective” length because it represents the number of possible alignment positions of a read to the transcript). That said, in practice it should make very little difference.

I see that in the post it is given EffectiveLength = Length - mean_fragment_length + 1; And how the mean_fragment_length is 300? One more question can I use this for Ensembl gene ids? or only for transcripts?

@beginner The fragment length is determined by your sequencing protocol. It’s (more or less) the distance between the two paired reads (including the read lengths, minus any overlap). For a typical long RNA-seq library prep it’s often in the order of 300–400 bp. You can run Picard’s CollectInsertSizeMetrics tool on your aligned BAM files to find the mean fragment length.

Sure. I will. I'm working with TCGA data. could you please tell about my other question I the previous comment.

You can use this on any annotated feature you like, both genes and transcripts. There are some fundamental issues with using genes (essentially, you lose information about alternative transcripts, and gene-level statistics by necessity need to average across transcripts) but that isn’t specific to quantification, and in practice it’s usually not an issue.

You can use countToFPKM package.
This package provides an easy to use function to convert the read count matrix into FPKM matrix; following the equation in

The fpkm() function requires three inputs to return FPKM as numeric matrix normalized by library size and feature length:

counts A numeric matrix of raw feature counts.
featureLength A    numeric vector with feature lengths that can be
obtained using
biomaRt.
meanFragmentLength A numeric vector with mean fragment lengths,
which can be calculate with Picard
using CollectInsertSizeMetrics.

See https://github.com/AAlhendi1707/countToFPKM


text:
Receptor-ligand database
How can I find ligands of receptors (or vice versa)?
For example, I want to find the receptors of IL-8 (CXCL8).

CellPhoneDB is a publicly available repository of curated receptors, ligands and their interactions.
The interactions can be searched or downloaded.

Have you tried STRING?
It's a publicly available protein-protein interaction database, so it won't be a complete picture.
Ingenuity Pathways Analysis (IPA) is really nice but is a paid service. Ask your institution if they have a license for one.

Can you include a link to IPA please.

https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/

CellTalkDB is another option for a searchable ligand-receptor database. According to the webpage it is curated from other databases and the literature. Details are available in the publication.


text:
featureCounts segmentation fault on Arch Linux
I am encountering a segmentation fault when attempting to run featureCounts from subread-1.6.3 on even small test data. I installed featureCounts from SourceForge using the precompiled Linux binary.
Most existing help threads on this topic (1, 2, 3) show some splash screen output before the error, but I get a segmentation fault before any error is shown.
$ featureCounts -T 1 -p -t exon -g gene_id -a /data/sai/annots/musmus_cdna_all_ensembl92.gtf -o test.txt test.bam test2.bam
Segmentation fault (core dumped)

I resolved this problem by installing from source; I assume the precompiled binary was not compatible with the Arch Linux kernel.

I encountered this with several tools on arch. It looks like most tools are compatible with Ubuntu (or Debian in general) or biolinux

Perhaps the authors should note the distribution of the binary build rather than denoting it as "linux"... or perhaps I shouldn't be so naive as to think such a thing could exist ;)


text:
about dotplot legend meaning

Here is the code I used:
DotPlot(object = cluster, genes.plot = c("Hcn4","Shox2","Nkx2-5","Nppa","Tnnt2","Gja5","Tbx18","Tbx3","Isl1","Scn5a","Cacna1d"),  plot.legend=TRUE)

I was confused about the legend with the Negative number shown for avg.exp.scale. Does it mean that the cluster expressed the genes at low level or that gene was down regulated? (like the Shox2 in the "0" compared with "2").

The expression values for each gene are scaled / standardized by subtracting the genes mean expression and dividing by its standard deviation. A value of -1 would imply it's one standard deviation below the mean expression for that gene.

Hi, can I ask you another question about the findmarker function in Seurat(https://satijalab.org/seurat/pbmc3k_tutorial.html) ? Can you tell me it's that possible that the avg_logFC shows negative number? Is it the same explanation as the legend question's answer?

The average log fold change is something like log(x) - log(y) where x and y are the average expression of two different groups.  A negative log fold change would mean the average expression of the second group (y) is higher than the first (x)

Hi, thank you for your explanation. It really makes sense to me. And if the code'cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
print(x = head(x = cluster1.markers, n = 5))' to find all markers of cluster1. Does it mean that x is cluster1 and y is the total cell?


text:
Mutation detection using Varscan2 on RNA sequencing for estimating tumour clones with pyclone or other package
I would like to analyze my RNAseq profiles from bulk tissue samples (Paired-End, 50M reads/sample, tumour-normal pairs) with varscan2 to detect mutations.
Then I would like to use those detected mutations to estimate presence of different tumour clones (likewise pyclone on WGS data). Therefore my questions:

Is it possible to estimate number of clones and clonal abundancy from RNA seq profiles from bulk tissue? and if yes: are 50M reads/sample enough coverage.

Could you recommend which package to use for this, for example is pyclone suitable for this?

Thank you in advance for the advice!

I guess this is not possible, as mutation frequencies (ratio of reads with and without mutation) used in this cluster approach depend highly on transcriptional regulation in the different clones, and does not only depend on the presence or absence of the mutation.