text:
Python scripting with ete3 to query NCBI's Taxonomy: "sqlite3 Warning (can only execute one statement at a time)"
I am using this script:
import csv
import time
import sys
from ete3 import NCBITaxa

ncbi = NCBITaxa()

def get_desired_ranks(taxid, desired_ranks):
    lineage = ncbi.get_lineage(taxid)
    names = ncbi.get_taxid_translator(lineage)
    lineage2ranks = ncbi.get_rank(names)
    ranks2lineage = dict((rank,taxid) for (taxid, rank) in lineage2ranks.items())
    return{'{}_id'.format(rank): ranks2lineage.get(rank, '<not present>') for rank in desired_ranks}

if __name__ == '__main__':
    file = open(sys.argv[1], "r")
    taxids = []
    contigs = []
    for line in file:
        line = line.split("\n")[0]
        taxids.append(line.split(",")[0])
        contigs.append(line.split(",")[1])

    desired_ranks = ['superkingdom', 'phylum']
    results = list()
    for taxid in taxids:
        results.append(list())
        results[-1].append(str(taxid))
        ranks = get_desired_ranks(taxid, desired_ranks)
        for key, rank in ranks.items():
            if rank != '<not present>':
                results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
            else:
                results[-1].append(rank)

    i = 0
    for result in results:
        print(contigs[i] + ','),
        print(','.join(result))
        i += 1

    file.close()

The script takes taxids from a file and fetches their respective lineages from a local copy of NCBI's Taxonomy database.  Strangely, this script works fine when I run it on small sets of taxids (~70, ~100), but most of my datasets are upwards of 280k taxids and these break the script.
I get this complete error:
Traceback (most recent call last):
  File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 31, in <module>
    ranks = get_desired_ranks(taxid, desired_ranks)
  File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 11, in get_desired_ranks
    lineage = ncbi.get_lineage(taxid)
  File "/data1/lstout/.local/lib/python2.7/site-packages/ete3/ncbi_taxonomy/ncbiquery.py", line 227, in get_lineage
    result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
sqlite3.Warning: You can only execute one statement at a time.

The first two files from the traceback are simply the script I referenced above, the third file is one of ete3's.  And as I stated, the script works fine with small datasets.
What I have tried:

Importing the time module and sleeping for a few milliseconds/hundredths of a second before/after my offending lines of code on lines 11 and 31.  No effect.
Went to line 227 in ete3's code...
result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])

and changed the "execute" function to "executescript" in order to be able to handle multiple queries at once (as that seems to be the problem).  This produced a new error and led to a rabbit hole of me changing minor things in their script trying to fudge this to work.  No result.  This is the complete offending function:
    def get_lineage(self, taxid):
    """Given a valid taxid number, return its corresponding lineage track as a
    hierarchically sorted list of parent taxids.
    """
    if not taxid:
        return None
    result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
    raw_track = result.fetchone()
    if not raw_track:
        #perhaps is an obsolete taxid
        _, merged_conversion = self._translate_merged([taxid])
        if taxid in merged_conversion:
            result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
            raw_track = result.fetchone()
        # if not raise error
        if not raw_track:
            #raw_track = ["1"]
            raise ValueError("%s taxid not found" %taxid)
        else:
            warnings.warn("taxid %s was translated into %s" %(taxid, merged_conversion[taxid]))

    track = list(map(int, raw_track[0].split(",")))
    return list(reversed(track))

What bothers me so much is that this works on small amounts of data!  I'm running these scripts from my school's high performance computer and have tried running on their head node and in an interactive moab scheduler.  Nothing has helped.

Please [edit] your question and include the script you are using. Both so we don't need to go back and forth between different tabs, and because the linked script may be edited. Also so we can see your modifications, understand what lines the error messages refer to etc.

Since you have 280k taxids i think querying is not the way to go; instead download a flat file from NCBI taxonomy with all the data and process that

@Chris_Rands I am using a local copy of taxonomy,  ncbi=ncbi.Taxa() gets the DB.

I overlooked the fact that (for my data) about 1:50,000 results returned from my blasting yielded a ";" within the hit denoting two taxids and did not account for these in my script.
How a given researcher will respond to their own results having ";"-hits will probably vary (you may want to do extra analysis, etc.).  In my case, I am only interested in the first taxid before the ";".  So I added a simply python .split() to my script before appending the taxid to the list that I query:
if __name__ == '__main__':
    file = open(sys.argv[1], "r")
    taxids = []
    contigs = []
    for line in file:
        line = line.split("\n")[0]
        tID = line.split(",")[0]
        if ';' in tID:
            tID = tID.split(";")[0]
        taxids.append(tID)
        contigs.append(line.split(",")[1])

    desired_ranks = ['superkingdom', 'phylum']
    results = list()
    for taxid in taxids:
        results.append(list())
        results[-1].append(str(taxid))
        ranks = get_desired_ranks(taxid, desired_ranks)
        for key, rank in ranks.items():
            if rank != '<not present>':
                results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
            else:
                results[-1].append(rank)

    i = 0
    for result in results:
        print(contigs[i] + ','),
        print(','.join(result))
        i += 1

    file.close()


text:
Creating a computational table to show the number of available genes between all possible pairs of combinations in eight gene lists
I am looking to create a table (on RStudio) that will show the overlap between 8 different gene lists wherever possible. I have saved the 8 lists as .txt files all in one folder and set that particular folder as my working directory. Any idea of code that will work for this?

If anyone knows what I can alter to make the current code work or suggest an alternative way I would be very happy to hear any suggestions!
data.file1 <- "GSE108363_BCGdown_V.txt"
data.file2 <- "GSE108363_BCGdown_D.txt"
data.file3 <- "GSE108363_BCGup_V.txt"
data.file4 <- "GSE108363_BCGup_D.txt"
data.file5 <- "GSE108363_MTBdown_V.txt"
data.file6 <- "GSE108363_MTBdown_D.txt"
data.file7 <- "GSE108363_MTBup_V.txt"
data.file8 <- "GSE108363_MTBup_D.txt"

genevect1 <- scan(data.file1, what=character(), sep="\n")
genevect2 <- scan(data.file2, what=character(), sep="\n")
genevect3 <- scan(data.file3, what=character(), sep="\n")
genevect4 <- scan(data.file4, what=character(), sep="\n")
genevect5 <- scan(data.file5, what=character(), sep="\n")
genevect6 <- scan(data.file6, what=character(), sep="\n")
genevect7 <- scan(data.file7, what=character(), sep="\n")
genevect8 <- scan(data.file8, what=character(), sep="\n")

filelist <- list(data.file1, data.file2, data.file3, data.file4, data.file5, data.file6, data.file7, data.file8)
#-------------------------------------------------------------------------------------------------------
# "Lapply" is basically making a "for loop".
#-------------------------------------------------------------------------------------------------------
all_gene_vectors <- lapply(filelist, scan, what=character(), sep="\n")

#-------------------------------------------------------------------------------------------------------
# Then make the intersection of *all* genesets
#-------------------------------------------------------------------------------------------------------
final_inter <- all_gene_vectors[[1]]
for (next_genevect in all_gene_vectors[2:length(all_gene_vectors)]) {
  final_inter <- intersect(final_inter, next_genevect)
}
show(final_inter)

#-------------------------------------------------------------------------------------------------------
# Or compute all pairwise intersections (untested)
#-------------------------------------------------------------------------------------------------------
pairwise_inters <- list()
#-------------------------------------------------------------------------------------------------------
# Iterate over all possible pairs i,j (i ≠ j)
#-------------------------------------------------------------------------------------------------------
N <- length(all_gene_vectors)
for (i in 1:(N-1)) {
  ith_inters <- list()
  for (j in (i+1):N) {
    ith_inters[[ names(all_gene_vectors)[j ]]] <- intersect(all_gene_vectors[[i]], all_gene_vectors[[j]])
  }
  pairwise_inters[[ names(all_gene_vectors)[i] ]] <- ith_inters
}

#-------------------------------------------------------------------------------------------------------
#Generate table in wide format:
#-------------------------------------------------------------------------------------------------------
#Concatenate gene names of all tables to a separate file
#Extract unique values in this file
#Perform outer join between the separate file and each table

all_names <- c(data.list)
for (tp in timepoints){
  all_names <- append(all_names, tables[[tp]]$genes)
}
all_names <- unique(all_names)

wide_table <- matrix(nrow=length(all_names), ncol=(1+length(timepoints)))
wide_table <- as.data.frame(wide_table)
colnames(wide_table) <- append("genes", timepoints)

for (tp in timepoints){
  tmp <- tables[[tp]]
  tmp$tp <- 1
  tmp <- merge(wide_table, tmp, by="genes", all=T)
  wide_table$tp <- tmp$tp
}

set.seed(11)
BCG_validation_Up <- sample(letters[1:429], 20)
BCG_discovery_Up <- sample(letters[1:250], 20)
MTB_validation_Up <- sample(letters[1:286], 20)
MTB_discovery_Up <- sample(letters[1:128], 20)
BCG_validation_Down <- sample(letters[1:267], 20)
BCG_discovery_Down  <- sample(letters[1:350], 20)
MTB_validation_Down <- sample(letters[1:244], 20)
MTB_discovery_Down <- sample(letters[1:86], 20)

cross_table <- matrix(, nrow = 4, ncol = 4)
rownames(cross_table) <- c("BCG_validation_Up", "BCG_discovery_Up", "MTB_validation_Up", "MTB_discovery_Up")
colnames(cross_table) <- c("BCG_validation_Down", "BCG_discovery_Down", "MTB_validation_Up", "MTB_discovery_Up")

for (i in 1:4){
  for(j in 1:4){
    cross_table[i,j] <- length(intersect(get(paste0("",i,"_Up")),(get(paste0("",j,"_Down")))))
  }
}

cross_table

Please show how the end table should look like, and give a small example to test with.

It looks like you already have code. Can you describe why what you have is not acceptable?

@winni2k Yes but I get the following error code when I attempt to execute it: Error in ith_inters[[names(all_gene_vectors)[j]]] <- intersect(all_gene_vectors[[i]],  :
attempt to select less than one element in OneIndex

@b.nota I have inserted a screenshot of the table I would like to be filled in, please find it above

That's a confusion matrix of 4 x 4, but you have 8 lists. Do you want a 8 x 8 confusion matrix or do you combine some of the lists?

Please don't post screenshots of code! That just makes the question heavier to load, useless for people using screen readers and, most importantly, makes the code you show useless to us. We need to be able to copy your code and run it locally to test any solutions we come up with. Please [edit] your question and include the code here directly. Use the formatting tools to format it as code.

@b.nota What i meant was a 4x4 table with 4 of them referring to upregulated genes and the other 4 for down regulated genes

@terdon Thanks for letting me know, I didn't know that it was so inconvenient for you and others. Please find my question edited.

@user3762, that was absolutely not clear from the question!

@b.nota Apologies, I hope it is all clear now

@b.nota Thanks for your suggested answer. I will check to see it works fine. Just wanted to ask why you have used (,20) in particular in all of the brackets? And why 11 in set.seed?

I made a workable example by generating lists of random letters (20 letters each). The set.seed command is to get the same letters each time. I made this workable example because you did not give a workable example.

Hi @user3762, I am sensing a bit of confusion and frustration in the comments to your question. I think some of that could be avoided if you edited your post to be about a Minimal, Complete, and Verifiable example as described here: https://stackoverflow.com/help/mcve

To make a confusion matrix of 4 lists against 4 other lists, you need to get the length of all intersecting genes for each pair. Here an example with letters instead of gene names.
set.seed(11)
list1_up <- sample(letters[1:26], 20)
list1_down <- sample(letters[1:26], 20)
list2_up <- sample(letters[1:26], 20)
list2_down <- sample(letters[1:26], 20)
list3_up <- sample(letters[1:26], 20)
list3_down <- sample(letters[1:26], 20)
list4_up <- sample(letters[1:26], 20)
list4_down <- sample(letters[1:26], 20)

cross_table <- matrix(, nrow = 4, ncol = 4)
rownames(cross_table) <- c("list1_up", "list2_up", "list3_up", "list4_up")
colnames(cross_table) <- c("list1_down", "list2_down", "list3_down", "list4_down")

for (i in 1:4){
  for(j in 1:4){
    cross_table[i,j] <- length(intersect(get(paste0("list",i,"_up")),(get(paste0("list",j,"_down")))))
    }
}

cross_table
         list1_down list2_down list3_down list4_down
list1_up         15         14         15         15
list2_up         18         15         15         15
list3_up         17         15         14         16
list4_up         17         15         15         16

Please could you look at my attempt at altering this code to suit my data and see where I have gone wrong? I have added the code to my question above

No, I will not look at your attempt. Here is my suggestion, do with it what you want.

@user3762 It would be better if you describe what has gone wrong and how the solution here fails you. Also show your attempts to fix or adapt this solution into one that works for you. Good luck!

@benn nice manners

source: Stack Exchange

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

sort: votes


text:
Why the t-test for a specific gene shows different value compared to differential analysis?
I have RNA-Seq data for LUNG cancer. 370 tumor and 50 Normal. For differential analysis initially I did some filtering and kept approx. 19k genes for further analysis. I used edgeR.
With a FC > or < 2 and FDR < 0.05
tr <- glmTreat(fit, contrast=contrast.matrix, lfc = log2(2))
I got approx 3000 differentially expressed genes.
Among the 3000 differentially expressed genes I'm interested in a gene OIT3.
GeneSymbol  logFC    unshrunk.logFC   logCPM          PValue      FDR
  OIT3  2.156487448  2.171599487    -0.437234244    0.000701037 0.005304312

I see the above gene OIT3 is Upregulated. Why I see logCPM with negative value?
And I plotted the logCPM expression of the above gene and did a t-test which gave pvalue 0.1 which means not significant.

1) I see the gene is differentially expressed with Differential analysis using edgeR. Why I see logCPM with negative value?
2) With t-test why it is not significant?

What is your design model and what is your contrast? It would help to clarify what is going on

If a t-test gave the same results as the more sophisticated software, why would anyone use the sophisticated stuff?  At first glance, logging your data is going to knock down those high outliers, and I bet edgeR is using a normalization of raw counts, not logged ones.

The $log(CPM)$ of any low-moderately expressed gene will be negative. There is nothing unexpected there.
Your statistics are inappropriate for a variety of reasons. Firstly, a CPM is not a robust value that's comparable between samples (this is why CPMs aren't used for statistics). edgeR performs more appropriate normalization and incorporates that into its model. Secondly, though you have a fair number of replicates, edgeR is pooling information from all of the genes to better estimate the variance.

I encourage you to read the edgeR paper and likely that of Limma too.

So, if I want to show boxplot of OIT3 gene. I should use only log2CPM on y-axis but should not calculate t-test with that. So, now In the above plot I should remove t-test pvalue and keep it is significant based on FC >2 and FDR < 0.05. Am I right?

If you want to show a boxplot you should use normalized counts.

Ok thanq very much for the answer. With edgeR how to get normalized counts?

There is no reason your t-test should reproduce edgeR. In fact, edgeR exists because t-test is inappropriate.
edgeR does the tests by pooling information from all genes, because with the low number of replicates your t-test doesn't give sufficient statistical power.
What you need to do is: check visually your gene and make a decision yourself based on graphical results and edgeR results if your gene shows biologically significant.

So, if I want to show boxplot of OIT3 gene. I should use only log2CPM on y-axis but should not calculate t-test with that. So, now In the above plot I should remove t-test pvalue and keep it is significant based on FC >2 and FDR < 0.05. Am I right?

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


text:
Can't install newest Blast from Conda
I have run into an odd situation, where Anaconda repo page clearly states that the newest version of Blast in Bioconda channel is 2.7.1, however the version Conda wants to install (command conda install -c bioconda blast) is 2.6.0.
Even specifying the version with command conda install -c bioconda blast=2.7.1 does not work and Conda returns a PackagesNotFoundError.
By observing the webpage mentioned above, one can see that there is a "v" before the version number, however adding this letter to the command with specified version did not help.
To make things even more odd, conda search blast returns a list of Blast versions available and there is a 2.7.1 version available in Bioconda channel.
I have added the Bioconda channel into my Conda.
My question is: How can I install the 2.7.1 (by the time of writing this question the newest version available in Anaconda repo) through Conda?
Edit: I have Ubuntu Server 18.04 LTS and Conda 4.5.11.

I was able to reproduce your problem, even in a blank environment.
$ conda create -n test
$ conda activate test
$ conda install bioconda::blast=2.7.1
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - bioconda::blast=2.7.1
  - boost=1.64
  - bioconda::blast=2.7.1
  - gnutls
  - bioconda::blast=2.7.1
  - nettle[version='3.3|3.3.*']
  - bioconda::blast=2.7.1

Note that boost, gnutls, and nettle are listed as requirements here.
If you do a search for those packages, you'll see that boost is available in a default channel, whereas gnutls and nettle are only in the conda-forge channel, which is not a default.
Adding those packages to the install command resolved this issue, for me.
$ conda install bioconda::blast=2.7.1 conda-forge::gnutls conda-forge::nettle --dry-run
Solving environment: done

## Package Plan ##

environment location: /Users/hawleyj/anaconda3/envs/test

added / updated specs:
    - bioconda::blast=2.7.1
    - conda-forge::gnutls
    - conda-forge::nettle
...

As mentioned in the comment below, to avoid this issue entirely, before using Bioconda packages make sure you run the proper channel configuration first.
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge

The appropriate channel order is mentioned in the documentation: https://bioconda.github.io/#set-up-channels


text:
How to create a graph of multi column gene dataset?
well, I'm new to Bio. what I'm going to do is create a network/graph of a dataset related to the genes I think. here's the sample of the nodes dataset with CSV format named lung-nodes :

and this is a sample of edges I think with trrust name :

these data are stored in two separate file
what I must do is create a graph(network) of these nodes and edges by Python and a visualization library like networkX, but I do not how.
I must make a relation between these to the dataset by creating a graph of them
can you help me out?

You must work in Python? Is this a homework assignment?

@benn, yes it is. I am not familiar so much with gene and protein properties. I already create a graph from social network dataset but they weren't like this. they had just two column in a single file. but right now i got two separate CSV file which contains bio info which i must create a graph of them. i think R studio could make it easier

As far as I understand, it seems that the first table you provided contains biological/experimental information related to your genes, i.e. the nodes of the network you will build. The second table is the one you need to build the network: each line is an edge. You can use it as input and build the network using the IGraph library which works for both Python (as you requested) or R (as you suggested as an alternative).
Most probably, but I can't say for sure, the first table was just a dataset from which you (or someone for you) recovered a list of nodes that where then submitted to a database, i.e. String, commonly used to recover information related to the interactions between your nodes, i.e. genes. Once you have this info (in the second table) you can build the network!
edit:
if you're not familiar with IGraph, you can have a look at this question describing how to obtain an IGraph graph from a csv file or this other one which uses a different csv file format to build a graph


text:
Subset data frame based on ID
I tried to subset the dataframe in R since there are empty space after each line in the sequence which shows 'NA' when i read the files into R ,as a result of which when i subset it just takes the sequence Protein_IDS and sequence line by line not as whole
Im looking for this Protein_ID in one column and respective sequence in other
Proteins_IDs                Sequence
  <chr>                       <chr>
1 fig|6666666.151743.peg.5284 MERQAVAVVFHDWRRAIQQRNKTAAADRGVDPVGVGHNSKIESLQPRSSVCRASLIQVPG
2 NA                          RLLNGEIFYSLAEAKVLIEAWRRHYNTVRPHSSLGYRPPAPETATPPYPASGSASLHLRP
3 NA                          DMAAMGLIH
4 NA                          NA
5 fig|6666666.151743.peg.5285 MIDRQTAQVRRWSFVLGSSIAAMAAASATPAWAQCSPDPTVADSVTNCTGTDSDGLRVTT
6 NA                          DGSMVEVAPSASVSGAGAAAITVDIAPVSGSYYARKATINVGGSVDGGTQSGISVLSGSV

So i have made another file with list of protein ID which i want to subset from this file..I tried removing empty space from the file but couldn;t as it the sequence should come in one line but instead there is line break
data files
I tried this SEQ <- inner_join(seq,protein_id) whose output is this its like each ID is associated with multiple sequence
fig|6666666.151743.peg.5044 MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAI
2 fig|6666666.151743.peg.5063 MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGT
3 fig|6666666.151743.peg.5135 MLRVADGYVEEAADDAIVPSPAVPLSVPSAPVRFEAAPARYGARRGPNVPAILVIVLLHA

It is very hard to understand your question. Please try to formulate the question, by telling what you have, what you want and what you have tried (including code).

edited the question

I want to advise you to use FASTA format, instead of obscure unknown formats. In R, there are many packages that can handle FASTA format and it is therefore more convenient. I will show you how to make a FASTA file from your data, and then I will show an example of how to import such a FASTA file and sub select your proteins of interest.
input <- scan("seq2.txt", what = character(), sep = "\t")
length(input)
new <- gsub("f", ">f", input[-1:-2])
write.table(new, "seq2.fasta", quote = F, row.names = F, col.names = F)

library(seqinr)
fasta.file <- read.fasta("seq2.fasta", seqtype = "AA", as.string = TRUE)
prot.id <- scan("protein_id.txt", what = character(), sep = "\t")
selection <- fasta.file[names(fasta.file) %in% prot.id]
length(selection)
[1] 87

head(selection, 2)
$`fig|6666666.151743.peg.5044`
[1] "MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAITSEQAIKLGITQPDVIAQMVPGLQLNHSAWGLTPFLRGVGTNSVAIGNEPAVALFVDDVYLPTGSAGVFQLNSISSIEVLKGPQGTLFGRNATGGVIHIRTRDPSFDPSASATASYGNYDTISGQLYGSTGLTDTVAINFALQGIDQRKGWGVNGTTGNEANSMSGWGGRAKLLWNMGEETSLLVAGSYDRRSSDQGQPGRSIAGTFTRDGYSADAAGLGFYDVMSNLDTVQKSTTRSLSGKLTHDFGSARLVLISAFSKTVVDPVVSDSDGSPSPFFQTREFAMGGRTFTQEVQLLSPEGSAVKWILGAFYLHDRSFMDATYEGQAFGTSKSHTDAVQRTNSIAGFGQVTTEILPRTNLTLGLRYTSDRRSLDAVARLLPATGAAVTFGPFEDAKTFNSVTGRASIDHHFSDDLMGYIAYNRGFKSGIYNFTGIVAGTTTLPPPVNPETLDAYTIGFKSEFLNRRLRVNAEAFLYEYSDLQVQIVLPTGTTLTNASKATIKGIDFEISAIPVRNLTVSASIGLLDGRYKDFRNAPSFFLLPPNAPTPYDPARCSGPNPYPAVQGSAPATQRSCDVSGNRLIQAPPFSSNLQLLYTLPTDLGDFDLSASWAHGGNYEYDVFSSPLTRQPRTDIVNGSVQWTAPGGTFNVRLWANNLTGEKYYTIATTGSTAGAKYAAAAPRTYGVTLGANF"
attr(,"name")
[1] "fig|6666666.151743.peg.5044"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5044"
attr(,"class")
[1] "SeqFastaAA"

$`fig|6666666.151743.peg.5063`
[1] "MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGTKREEPLARLPISVAIVQAPGEQGGALLPGTEEVAQRLDGLTLTNLGPGRNRQFIRGVADSPFNGNSQSTVAILINDARATYNAPDPDLKLVDVERVELLKGPQGSLYGTGTLGGVYHIVTRPPDAGSFAAMAGMGASLATHGEMGGSAEAMVNLPVVADRLALRAVAYAVREPGWIDTAGRGDANRGELRGGRVALRALAGDWTVDLNGMAQMRDVADSQYVYARHALERPAQPAEPHDNDFALMALRVAGKLGGADLVYAGNWVDHDVSSTLDASDLSASAFGIAPPALYLEQRRNRLFDQEVRLSQTRDDGLSWLAGLSFVDVRNHLSGVVRGVRDETVGTGIQGVVEFSAFGDMRMPLTGALSLSVGGRLFHSRIEEEQGARHARIEASVARSGFTPSLALSWQAGAHDFLYARYAGALRPGGLSRTGAASAERFESDELQTVEMGWRHSASAWLTVNLDAFLTTWSHMQSDYLLADGLIGTRNVGDARIRGVEASVDWKPAAGWAISAGLVGQSALLTRTNEGLRLEDRHLPVVPDWKARFGVNRHFDLSGWSGSMGVNLAYLGPSRLSFDPGLDRRMGDYALANADLALRRKDWSIRLSGSNLLDCKADSFAFGNPFALARGTQYTPLRPRQWAIGFTKNW"
attr(,"name")
[1] "fig|6666666.151743.peg.5063"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5063"
attr(,"class")
[1] "SeqFastaAA"

looks like its bit more complicated than i thought let me run ,do you have a sed or awk ways to get around this file format?

No, like I said: I advice you to USE the fasta file format because you can use a lot of tools for your purposes instead of reinventing the wheel with your own format.

how did you write the selection object into a file as the object is list .."Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :
cannot coerce class ‘"SeqFastaAA"’ to a data.frame"

able to write it...and im really glad for your help.....

Great to hear, good luck with the analysis.


text:
multiple COSMIC id for the same mutation
I would like to know why there are different entries for the same mutation in COSMIC, for example COSM6954941 and COSM12833 both refer to ERBB4 c.908C>A.
In this specific case the field Gene name in the two entries is different, COSM6954941 report ERBB4_ENST00000402597 while COSM12833 report ERBB4, but there are many cases in which the Gene name is the same, for example COSM1361 and COSM1219.

Welcome to the site! Please [edit] your question and add links to these cosmic entries so we can check what you mean more easily. You might not be able to link to all of them (there are some restrictions on how many links new users can include, to limit spam). If you have issues with that, just add the links as text (remove the "https://:") and I will add them for you.

Hi @terdon, I added the links as you requested. Please note that I ave got an answer via email that I report using the "answer your question" button.

Great! And thank you for remembering to come back and post the answer. That's very helpful and might help future users with the same problem.

Harry Boutselakis addressed the question via the customer support of the COSMIC web site.

Dear Ivan, Thank you for using COSMIC.

I would like to know why there are different entries for the same mutation   in COSMIC, for example COSM6954941 and COSM12833 both refer
    to ERBB4   c.908C>A.
In this specific case the field Gene name in the two entries is different,   COSM6954941 report ERBB4_ENST00000402597 while COSM12833
    report ERBB4,

The mutations are annotated onto two different isoforms of ERBB4. So
  this annotation refers to two different transcripts. In these cases
  the consequence of the mutation on the mRNA and polypeptide is the
  same reflected by the syntax.

but there are many cases in which the Gene name is the same, for example   COSM1361 and COSM1219.

This case the consequence one mutation on the mRNA is unknown , but
  the amino acid change is the same on both cases.
The information stored in COSMIC is curated from experimental data,
  depending on the evidence and experimental methodology used, the
  mutations are differentiated accordingly. I hope this answers your
  question.
Kind regards,
Harry


text:
bcftools filtering all files in a directory
Probably a silly oversight on my part, but I'm trying to filter all the vcfs in a directory with bcftools using a simple loop. My basic command is working fine:
bcftools filter -i 'QUAL > 1000' -o filter/file1out.vcf file1in.vcf

but when I try to loop it, it echoes through each filename, but generates no output files or error messages.
for f in *.vcf; do echo "filtering $f"; bcftools filter -i 'QUAL > 1000' -o filter/$f_out.vcf $f; done

Am I not handling the output direction properly? Any help much appreciated!

I assume your shell is looking for the variable $f_out, but can't find it because you did not define it. It probably generated a hidden file with the name .vcf in your output directory.
Instead, use the following: -o ${f}_out.vcf. This will "protect/limit" the bash variable to just $f rather than $f_out.
When encountering issues like this it's worth putting echo in your loop right before your command, in this case bcftools (adding quotes whenever necessary). That would show you the command as how it is going to be executed and would show you, in this case, that the output name is incorrect.

This works thanks! But I wonder why I've had a similar loop working fine with $f.SnpSift.vcf as my output? Does it understand the . differently than the _ after $f?

Yes that's right, and easy to test using a dummy variable and echo

Alternatively, instead of (or in addition to) echo-ing you can use set -u which will give an error on undefined variables (this assumes you use Bash or Zsh as your shell (not sure about others)).

Using a for loop for this is by far the safer option. However, you need to quote your variables and also make sure you're passing the right variable name. While you're at it, you can also remove the extension so that your filtered file name can be foo_out.vcf instead of foo.vcf_out.vcf (which is what your loop would create. Putting all that together:
for f in *.vcf; do
    echo "filtering $f";
bcftools filter -i 'QUAL > 1000' -o filter/"${f//.vcf}"_out.vcf "$f"; done
fi

I like the "${f//.vcf}", thanks - my filenames have been growing unnecessarily long!

For one-liners, xargs is often more convenient than the for loop:
ls *.vcf | sed s,.vcf$,, \
  | xargs -i echo 'bcftools filter -i "QUAL > 1000" -o "filter/{}_out.vcf" "{}.vcf"' \
  | parallel

Note that your command line (with @WouterDeCoster's fix) writes output to foo.vcf_out.vcf. My command line writes to foo_out.vcf.

Xargs is indeed convenient, but parsing ls is very fragile and will break on any non-standard file name. Why not for f in *vcf; do printf 'bcftools filter -i "QUAL > 1000" -o filter/"%s" "%s"\n' "${f//.vcf}_out.vcf" "$f"; done | parallel? That seems as simple as xargs but can deal with arbitrary file names.

@terdon I don't see that is an issue. You can simply add quotation marks.

Yes, that will protect from spaces, but not from newlines. And OK, having newlines in filenames isn't something that happens very often but it is possible. I freely admit, however, that I may just be nitpicking because I'm too used to hanging out on [unix.se] where this sort of thing is considered very important and any example of parsing ls is met with horror. Still, since the shell loop is safe and painless, it just seems better to use that instead.


text:
Question about research in bioinformatics for Computer Science student
I am a computer science student, who was interested in taking research in the field of bioinformatics.
The first research idea I had was to make a model that could predict the possibility of type 2 diabetes by using machine learning such as convolutional neural networks or recurrent neural networks by training the model with fasta data of protein that causes diabetes.
The research will focus more on the architecture of the machine learning itself.
Does this research possible to work on ? If yes, how and where can I get the data?
All my lecturers have no background in bioinformatic, so I can't ask any of them.
Edit:
I m sorry, not protein but gene like CAPN10, CDKAL1, CDKN2A, HHEX, HNF4A, IGF2BP2, KCNJ11, PPARG, SLC30A8, and TCF7L2 that shown have risk to be associated in Type 2 Diabetes (http://diabetes.diabetesjournals.org/content/56/12/3033).

Before you can do anything valuable with machine learning in bioinformatics you must understand the biology. You cannot select meaningful features if you do not know what they mean. There is no protein that causes diabetes, and if there was its fasta sequence would most likely not be a very helpful feature.

I m sorry, not protein but gene like CAPN10, CDKAL1, CDKN2A, HHEX, HNF4A, IGF2BP2, KCNJ11, PPARG, SLC30A8, and TCF7L2 that shown have risk to be associated in Type 2 Diabetes (http://diabetes.diabetesjournals.org/content/56/12/3033). The model will be trained by those gene, both from the infected and healthy sample,it is still possible? Thank you for answering

Are you familiar with GWAS, linkage disequilibrium, functional vs meaningless variants, sequencing and variant calling? If not, that would be your starting point. If so, you would know that the variants identified are most likely non-coding and not the true functional variant. Follow-up studies, ideally including functional modeling are necessary to figure out those variants which actually affect diabetes and which are just on the same haplotype.

In addition, while I don't work in diabetes research I would expect you can find a GWAS with a larger number of individuals than the one you linked. Ideally, you would use a dataset from tens of thousands or more individuals.

I think , this kind of research never stop until find way for cure .
I have some research in heart disease,and write this paperECG, you can have some research in improve learning method to find fast and more accurate diabetes , there are lot of database your can search in google and more famous is UCI   .


text:
What does an FDR value of 1 in RNA-seq mean?
I am looking at the supplemental data from the paper "An allelic series of miR-17 ∼ 92-mutant mice uncovers functional specialization and cooperation among members of a microRNA polycistron" which lists the genes that are differentially expressed between a particular knock-out mouse and the wild-type control. In addition to showing the logFC changes for each of the genes, the table includes their FDR value. In many cases, the FDR value is 1. What does this mean?

A link/doi to the paper would be helpful.

Sure! Here is the link: https://www.ncbi.nlm.nih.gov/pubmed/26029871 Thank you very much, Leah

@leah I don't see where they are reporting FDR values of 1 in that paper, but Mowgli's answer is correct for what FDR is and likely correct that the numbers you are seeing are meant to be FDR corrected p-values. They shouldn't be called "FDR values" though so the authors might have made a mistake or just an unfortunate choice of abbreviation (or you did). FDR is a threshold, sort of like when you choose an alpha value as a criterion for significance.

FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing.
As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there's a 5% chance that what you call a hit is actually not signal, just some extreme value that happened by chance, and really there's nothing to see here".
This definition has deep consequences if you run many experiments: imagine that you run 100 experiments (say, you test the effect of 100 different molecules), if none of them has any effect, on average you could still get 5% of false positives. So if you run these 100 experiments, and you get, say, 10 hits with a p-value below your thresold of 0.05; how can you tell which ones are real and which ones are just chance in the absence of biological signal?
This is where you can use FDR.
If you control for a false discovery rate of, say, 0.2 (20%), it means that after a few computations (the Benjamini-Hochberg procedure described on wikipedia is simple and commonly used) you will lower your p-value cutoff so as to be sure that, among the N hits (probably less than the 10 you initially selected above) that you got above, a maximum of 20% are actually false positives (i.e. no biological effect, but surprisingly high/low signal values by pure chance).
In the article you mention they have a different FDR value for each sample, which is unusual. The values are not round numbers, so it seems likely that they are numbers obtained after applying the FDR correction (i.e. they are "transformed p-values") on which they based their decision to include the samples or not. The FDR cutoff that they chose should be written somewhere in the main text.

In practice, you can interpret it just like a p-value.
An FDR value is a p-value adjusted for multiple tests (by the Benjamini-Hochberg procedure). It stands for the “false discovery rate” it corrects for multiple testing by giving the proportion of tests above threshold alpha that will be false positives (i.e., detected when the null hypothesis is true). Note this is less stringent than the Holm-Bonferroni FWER adjusted p-value (which considers the chance of any false positives among significant results).
So for a FDR p-value of 0.05, up to 5% of these tests will be false positives (note this is not the case for uncorrected raw p-values).
For an FDR p-value of 1, up to 100% of these tests will be false positives. This makes sense since if you take test with a p-value if 1, then it should include all of the negative results, where the null hypothesis cannot be rejected. Only tests with p-values lower than this, have any chance of being rejected under the null hypothesis as true positives. If you are getting FDR adjusted p-values, you have done so many tests that you do not have the power to detect true positives and cannot rule out these being false positives due to noise of multiple testing. So an FDR p-value of 1 is definitely not significant under any circumstances.

I think you cannot interpret it exactly like a p-value in a Benjamini Hochberg procedure. It requires 1) ranking the p-values, 2) computing adjusted values based on rank and 3) finding the largest index satisfying the FDR inequality computation and 4) selecting all the values below this index. There may be values that end up being selected (ie. their index is sufficiently small) even though the computation leads to a value above your selection threshold.

Additionally, an FDR just tells you about the maximum percentage of false positives that you are willing to accept - so even correcting for an FDR of 1 can lead to 100% of correct hits, it just means that you are willing to accept up to 100% of type I errors in a worst-case scenario (and you just don't know if it happens or not)

For the case of p = 1 or in terms of interpreting the results of an experiment, it’s not much different. Any algorithm to compute FDR correction will compute them in order for you. It’s not necessary for a biologist to know the details, although this may be more appropriate for Cross Validated.

I don't quite agree; this nuance can add tens of hits to your results in a large scale experiment like high-throughput RNAseq. And yes, I think it's necessary for a biologist - or any other scientist, for that matter - to understand precisely the statistical tools they use for research :) my personal opinion.

For this specific case, it is a question about FDR of 1. This should not be used as a threshold for defining hits under any circumstances. This will make every result a hit by definition. It’s the most false positives you will allow and you are allowing all them. Researchers should understand the tools they’re using but they should seek further information on them, this is not the right forum for that level of detail.


text:
De novo motif discovery in protein sequences
I am trying to build a features matrix to be used for Random Forest based classification. I'd like to add, as features, short motifs which are common to all the protein sequences belonging to a specific gene family.
I tried to use MEME but I'd like to know if there are also other good tools (online or command line). Keep in mind I have more than 10.000 sequences and that I prefer short motifs (3/4 aa long).

I think MEME is a good tool for your purpose, but there are others as well. I can think of InterProScan for example (although I am not sure if it is really de novo). Here a summary of available tools.


text:
How to find genomic overlaps faster?
I was trying to find a way to find overlaps between two genomic ranges. I found a post at Biostars but it couldn't consider the chromosome information. For example:
library(IRanges)
df1 = data.frame(chr = c("chr1", "chr12"), start = c(10000, 10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

It returns the overlap from both chromosomes. I am wondering if there is a way to keep the chromosome information intact.

Just wanted to follow up to a question asked at Biostars. My solution and suggestion when dealing with genomic ranges.
library(GenomicRanges)
df1 = data.frame(chr = c("chr1",  "chr12"), start = c(10000,  10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

# define a genomic interval tree (git) for faster search

df1.gr = GRanges (IRanges(start = df1$start, end = df1$end), seqnames=df1$chr)
df2.gr = GRanges(IRanges(start=df2$posn, end = df2$posn), seqnames = df2$chr)

df1.git  = GNCList(df1.gr)
t1 = Sys.time ()
overlap_git = findOverlaps(df2.gr, df1.git)
t2 = Sys.time()
sub1 = difftime(t2, t1, tz, units = c("auto"))
print (sub1)
print (overlap_git)

t1 = Sys.time ()
overlap_intersect = intersect(df1.gr, df2.gr)
t2 = Sys.time ()
print(overlap_intersect)
sub2 = difftime(t2, t1, tz, units = c("auto"))
print(sub2)

One should be informed about differences between IRanges and GRanges (for example chromosome context); for example, findOverlaps with IRanges returns intersections from both chr1 and chr12, but with GRanges it doesn't
Interval Tree based search is way faster (almost 10x) than simple intersect function


text:
Linking of MeSH terms and genes
Is there an alternative to mesh2gene for linking MeSH terms to genes?
The http://www.ncibi.org/gene2mesh.html site appears to be down for months now (click on "Launch Gene2MeSH via gene2mesh.ncibi.org") and I've been unable to find an alternative so far...

I think you could use the MESH.db package in Bioconductor to annotate your genes. I have never used, but you should be able to do something like:
library("MeSH.db")
select(MeSH.db, keys = myGenes, keytype = "ENTREZ", columns = "MESH")

It depends on the MeSHDbi package from which you could retrieve other data form that database if needed


text:
Bash script error at paste command
I wrote script for pasting rsids on CADD output. Here is script.
#!/bin/bash
cd tmp
cut -f 1,2 CADD.tsv > fileA
paste fileA <(cut -f 2,125  CADD.tsv) > myNewFile
bedtools intersect -a  myNewFile -b New.vcf -wb |cut -f 1-4,7 > CADD.rsids.tsv

I have tested commands one by one and they work fine. But when I run script of these commands, it gives me this error.
$ sh cadd.rsids.sh
4: cadd.rsids.sh: Syntax error: "(" unexpected

Can anyone please tell me how to resolve this error?
Here are few lines of CADD.tsv file.
"##CADD GRCh38-v1.4 (c) University of Washington, Hudson-Alpha Institute for Biotechnology and Berlin Institute of Health 2013-2018. All rights reserved."
#Chrom  Pos Ref Alt Type    Length  AnnoType    Consequence ConsScore   ConsDetail  GC  CpG motifECount motifEName  motifEHIPos motifEScoreChng oAA nAA GeneID  FeatureID   GeneName    CCDS    Intron  Exon    cDNApos relcDNApos  CDSpos  relCDSpos   protPos relProtPos  Domain  Dst2Splice  Dst2SplType minDistTSS  minDistTSE  SIFTcat SIFTval PolyPhenCat PolyPhenVal priPhCons   mamPhCons   verPhCons   priPhyloP   mamPhyloP   verPhyloP   bStatistic  targetScan  mirSVR-Score    mirSVR-E    mirSVR-Aln  cHmm_E1 cHmm_E2 cHmm_E3 cHmm_E4 cHmm_E5 cHmm_E6 cHmm_E7 cHmm_E8 cHmm_E9 cHmm_E10    cHmm_E11    cHmm_E12    cHmm_E13    cHmm_E14    cHmm_E15    cHmm_E16    cHmm_E17    cHmm_E18    cHmm_E19    cHmm_E20    cHmm_E21    cHmm_E22    cHmm_E23    cHmm_E24    cHmm_E25    GerpRS  GerpRSpval  GerpN   GerpS   tOverlapMotifs  motifDist   EncodeH3K4me1-sum   EncodeH3K4me1-max   EncodeH3K4me2-sum   EncodeH3K4me2-max   EncodeH3K4me3-sum   EncodeH3K4me3-max   EncodeH3K9ac-sum    EncodeH3K9ac-max    EncodeH3K9me3-sum   EncodeH3K9me3-max   EncodeH3K27ac-sum   EncodeH3K27ac-max   EncodeH3K27me3-sum  EncodeH3K27me3-max  EncodeH3K36me3-sum  EncodeH3K36me3-max  EncodeH3K79me2-sum  EncodeH3K79me2-max  EncodeH4K20me1-sum  EncodeH4K20me1-max  EncodeH2AFZ-sum EncodeH2AFZ-max EncodeDNase-sum EncodeDNase-max EncodetotalRNA-sum  EncodetotalRNA-max  Grantham    Dist2Mutation   Freq100bp   Rare100bp   Sngl100bp   Freq1000bp  Rare1000bp  Sngl1000bp  Freq10000bp Rare10000bp Sngl10000bp EnsembleRegulatoryFeature   dbscSNV-ada_score   dbscSNV-rf_score    RemapOverlapTF  RemapOverlapCL  RawScore    PHRED
1   3362704 T   A   SNV 0   Transcript  INTRONIC    2   intron  0.536423841 0.04    NA  NA  NA  NA  NA  NA  ENSG00000142611 ENST00000270722 PRDM16  CCDS41236.2 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  118572  24215   NA  NA  NA  NA  0   0   0   -0.405  -1.689  -1.593  964 NA  NA  NA  NA  0   1   1   5   0   0   1   0   0   0   0   1   0   1   12  0   0   0   0   0   8   11  5   2   0   NA  NA  2.46    -4.93   1   0.02    4.53362 1.5498  3.49082 1.00224 9.30658 3.84201 5.88121 2.95189 9.31556 5.27366 4.01294 1.08327 25.99   5.10373 4.47027 1.51567 1.55624 0.7046  4.95028 1.49632 8.57929 1.45764 1.19689 0.644106    0.01783 0.01476 NA  16  1   0   15  4   10  184 17  101 1741    NA  NA  NA  5   19  -0.173019   0.424
1   7785635 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049245 ENST00000054666 VAMP3   CCDS88.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7785635 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7803233 C   T   SNV 0   Transcript  INTRONIC    2   intron  0.331125828 0.013333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  182 726 NA  NA  NA  NA  0.014   0   0   -0.553  -1.498  -0.938  673 NA  NA  NA  NA  0   0   0   6   0   0   23  5   0   0   0   0   0   0   14  0   0   0   0   0   0   0   0   0   0   NA  NA  7.49    -15 NA  NA  4.36329 1.02663 4.55863 1.4118  4.72226 1.51608 6.80823 2.09997 7.44497 2.86816 5.71457 2.35876 3.06623 0.62198 18.3251 3.23263 2.1423  0.84413 8.63178 2.70111 6.29127 1.59229 0.198671    0.048037    0.0465  0.01705 NA  6   1   1   15  2   2   134 26  39  1314    NA  NA  NA  NA  NA  -0.395046   0.067
1   7808665 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7808665 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7827519 C   G   SNV 0   CodingTranscript    NON_SYNONYMOUS  7   missense    0.529801325 0.066666667 NA  NA  NA  NA  P   A   ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 NA  18/22   2854    0.45172523  2590    0.712909441 864 0.714049587 lcompl  NA  NA  177 17277   tolerated   0.05    benign  0.127   0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 27  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   RegulatoryFeature   REGULATORY  4   regulatory  0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  NA  ENSR00000000832 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7828155 A   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.364238411 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  813 17022   NA  NA  NA  NA  0.009   0   0   -0.389  -0.735  -0.731  522 NA  NA  NA  NA  0   0   0   4   0   0   20  4   1   0   0   0   0   0   14  1   0   0   3   0   0   0   0   1   0   NA  NA  4.56    4.56    1   0.12    10.0753 2.35506 38.5512 20.107  21.923  10.3608 21.217  11.8934 5.61317 1.38406 9.74608 6.16783 3.76905 1.10216 24.6436 5.18488 14.8421 7.93141 4.56698 1.13989 18.7058 4.31099 0.500167    0.211141    0.16854 0.06818 NA  20  2   0   14  4   4   159 34  57  1421    NA  NA  NA  11  14  0.225659    4.402

Here are lines of New.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
1   3362704 rs11807862  T   A   923.01  PASS    "BaseCounts=31,0,0,30;BaseQRankSum=-0.108;DB;Dels=0;FS=3.561;GC=61.35;HaplotypeScore=1.7256;MQ=60;MQ0=0;MQRankSum=-0.368;QD=15.13;ReadPosRankSum=0.498;DP=282;AF=0.5;MLEAC=1;MLEAF=0.5;AN=14;AC=7"
1   7785635 rs228729    T   C   2294.01 PASS    "BaseCounts=0,44,0,0;DB;Dels=0;FS=0;GC=42.39;HaplotypeScore=0;MQ=60;MQ0=0;QD=28.13;BaseQRankSum=2.114;MQRankSum=-1.268;ReadPosRankSum=-0.548;DP=1430;AF=0.5;MLEAC=1;MLEAF=0.5;AN=70;AC=51"
1   7803233 rs228642    C   T   2082.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=-0.963;DB;Dels=0;FS=2.884;GC=44.89;HaplotypeScore=0.9999;MQ=60;MQ0=0;MQRankSum=0.058;QD=14.92;ReadPosRankSum=0.409;DP=1468;AF=0.5;MLEAC=1;MLEAF=0.5;AN=62;AC=44"
1   7808665 rs228666    T   C   1925.01 PASS    "BaseCounts=0,36,0,26;BaseQRankSum=0.064;DB;Dels=0;FS=0.979;GC=30.42;HaplotypeScore=1.7333;MQ=60;MQ0=0;MQRankSum=-0.82;QD=18.03;ReadPosRankSum=-1.32;DP=939;AF=0.5;MLEAC=1;MLEAF=0.5;AN=42;AC=25"
1   7827519 rs228697    C   G   786.01  PASS    "BaseCounts=0,14,24,0;BaseQRankSum=0.016;DB;Dels=0;FS=4.925;GC=58.35;HaplotypeScore=0;MQ=60;MQ0=0;MQRankSum=-0.698;QD=18.63;ReadPosRankSum=0.99;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=9"
1   7828155 rs2859388   A   G   2385.01 PASS    "BaseCounts=0,0,56,0;DB;Dels=0;FS=0;GC=35.41;HaplotypeScore=0;MQ=60;MQ0=0;QD=24.91;BaseQRankSum=-0.231;MQRankSum=0.948;ReadPosRankSum=-0.128;DP=1136;AF=0.5;MLEAC=1;MLEAF=0.5;AN=46;AC=31"
1   7830057 rs2640909   T   C   1421.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=1.051;DB;Dels=0;FS=5.041;GC=56.11;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=0.496;QD=16.61;ReadPosRankSum=-0.526;DP=638;AF=0.5;MLEAC=1;MLEAF=0.5;AN=32;AC=19"
1   11107089    rs12139042  G   A   1668.01 PASS    "BaseCounts=43,0,0,0;DB;Dels=0;FS=0;GC=45.64;HaplotypeScore=0.734;MQ=60;MQ0=0;QD=29.7;BaseQRankSum=-1.652;MQRankSum=0.027;ReadPosRankSum=-1.192;DP=354;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=10"
1   11787392    rs3737967   G   A   719.01  PASS    "BaseCounts=23,0,23,0;BaseQRankSum=-0.475;DB;Dels=0;FS=4.566;GC=61.85;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=-0.84;QD=15.63;ReadPosRankSum=0.329;DP=378;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11790870    rs2274976   C   T   777.01  PASS    "BaseCounts=0,18,0,24;BaseQRankSum=-0.042;DB;Dels=0;FS=4.677;GC=56.11;HaplotypeScore=0.6651;MQ=60;MQ0=0;MQRankSum=1.204;QD=15.05;ReadPosRankSum=0.836;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11792243    rs1476413   C   T   1751.01 PASS    "BaseCounts=1,18,0,21;BaseQRankSum=1.028;DB;Dels=0;FS=0;GC=56.11;HaplotypeScore=1.7287;MQ=60;MQ0=0;MQRankSum=-0.408;QD=16.65;ReadPosRankSum=-0.099;DP=954;AF=0.5;MLEAC=1;MLEAF=0.5;AN=50;AC=29"
1   11794400    rs4846051   G   A   2616.01 PASS    "BaseCounts=46,0,0,0;DB;Dels=0;FS=0;GC=59.35;HaplotypeScore=6.2573;MQ=60;MQ0=0;QD=31.75;DP=1654;AF=1;MLEAC=2;MLEAF=1;AN=78;AC=78"
1   11794419    rs1801131   T   G   2028.01 PASS    "BaseCounts=0,1,21,21;BaseQRankSum=-0.117;DB;Dels=0;FS=9.995;GC=59.35;HaplotypeScore=2.4328;MQ=60;MQ0=0;MQRankSum=1.03;QD=14.84;ReadPosRankSum=0.665;DP=1176;AF=0.5;MLEAC=1;MLEAF=0.5;AN=54;AC=33"

The error complains about the parenthesis "("

If I remove "(" parathesis, this gives me error of 4: c.rsids.sh: cannot open cut: No such file

Try to write a pipe with cut followed by paste command, see manual for paste command here.

I have used this command cut -f 2,125 CADD.tsv | paste fileA but it gives me only first two columns which are 1 and 2.

cut -f 2,125 CADD.tsv | paste fileA - I think you need to put a dash at the end?

Sorry for bothering again but putting dash at the end paste only column #2, it does not give column 125

To be honest, both lines work on my machine: paste fileA <(cut -f 2,125  CADD.tsv) > myNewFile and cut -f 2,125 CADD.tsv | paste fileA - they give me 4 columns, last one called PHRED.

My first guess is that the script isn't being run by bash. What operating system are you using? What exact command do you use to launch the script? Also, both of your files have headers which are not the same so should be skipped when pasteing.

I am using Ubuntu 16.04.5. I used this command sh Cadd.rsids.sh. I do not want header just columns. When I run these commands separately, they run without error. But not in script form. :(

Your problem is that you have written a script in bash but are then running it using sh. Bash, the Bourne-again shell, is not the same as the Bourne shell (sh). What is more, on Ubuntu, sh is actually not even the Bourne shell but another minimal shell called dash. The /bin/sh on Ubuntu is a symlink to dash:
$ ls -l /bin/sh
lrwxrwxrwx 1 root root 4 Mar  1  2018 /bin/sh -> dash

Dash is a POSIX compliant shell which means it has the features defined by the POSIX specification for the sh shell. In your case, the problem is caused by the fact that the <() input redirection construct is supported by bash but not by dash. So, in order to run your script, you can do one of two things:

Make the script executable (chmod a+x cadd.rsids.sh), and then run it directly:
/path/to/cadd.rsids.sh

This is the simplest approach since you already have the shebang line (#!/bin/bash) which will cause the script to be run by bash and not sh.
Run the script with bash instead of sh:
bash /path/to/cadd.rsids.sh

Instead of
sh /path/to/cadd.rsids.sh


text:
How to output all sequences with bwa mem, not `*`?
I've been running bwa mem -a for alignment, using the -a flag---this will

output all alignments for SE or unpaired PE

I've noticed in the SAM that there are several alignments with * in the SEQ and QUAL fields. Based on the documentation:

SEQ: segment SEQuence. This field can be a ‘*’ when the sequence is not stored. If not a ‘*’, the length of the sequence must equal the
  sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the
  base is identical to the reference base. No assumptions can be made on
  the letter cases.
QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A base quality is the phred-scaled base
  error probability which equals −10 log10 Pr{base is wrong}. This field
  can be a ‘*’ when quality is not stored. If not a ‘*’, SEQ must not be
  a ‘*’ and the length of the quality string ought to equal the length
  of SEQ.

it appears the sequence isn't stored.
I would strongly prefer to have the sequence in this case. Is there any to direct bwa to output these sequences?

Are you aware that many of these records are secondary/supplementary alignments? Does your workfow/model require/benefit from these alignments and handle them correctly?

@DanielStandage "Are you aware that many of these records are secondary/supplementary alignments?" Yes. " Does your workfow/model require/benefit from these alignments and handle them correctly?". Benefit, yes.

I don't believe that is an option.  Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.

Thanks. Did @user172818 ever give a reason why this feature isn't included?

see also https://www.biostars.org/p/352930/

@Pierre This is a great tool Pierre! I think it runs successfully for me. The problem is, I get errors downstream when using samtools view, i.e.. [E::sam_parse1] CIGAR and query sequence are of different length Are there work arounds?

I didn't test it with all the possible cases. Can you please try to isolate a minimal sam file containing the faulty reads ( all the reads with the same name) and send a bug report ? https://github.com/lindenb/jvarkit/issues/

@Pierre I may send you an e-mail instead. Thanks Pierre!

plindenbaum reverse(oohay) fr but I'll need a minimal failing SAM file.


text:
Real time transcript profiles
Do any methods exist (or are in the process of development) for investigating transcript data without lysing the cells, i.e, destroying the sample?

You may find this review useful: https://www.cell.com/molecular-cell/pdf/S1097-2765(18)30591-4.pdf

I'm not sure this question is about bioinformatics, because there isn't any method that is able to produce the data the way you ask. Or if it would exists it doesn't ask about any tool or method analysing/processing it

Answer:
There are no methods that exist for investigating transcript data inside a cell without destroying the cell. Separating RNA from cells requires essentially destroying the cells and using a chemical cocktail to split its constituent components.
Alternative 1: One option is to investigate extracellular transcripts of living cells is to centrifuge cells in culture and isolating RNA from the growth media supernatant [1]. Cell-free DNA and cell-free RNA are hot topics in contemporary biology and bioinformatics as a potential source for disease biomarkers. [2]
Alternative 2: Another option to monitor a transcript from a living cell is to use a transgenic construct with a reporter molecule like luciferase or a fluorescent protein to produce light or become fluorescent when a promotor is activated [3, 4]. This does not directly measure transcript quantity but answers a similar question: how does transcription change in a living cell?
Sources:
[1] Théry, C., Amigorena, S., Raposo, G., & Clayton, A. (2006). Isolation and characterization of exosomes from cell culture supernatants and biological fluids. Current protocols in cell biology, 30(1), 3-22.
[2] Yu, X. M., Wu, Y. C., Liu, X., Huang, X. C., Hou, X. X., Wang, J. L., ... & Ling, Z. Q. (2016). Cell-Free RNA content in peripheral blood as potential biomarkers for detecting circulating tumor cells in non-small cell lung carcinoma. International journal of molecular sciences, 17(11), 1845.
[3] Kaskova, Z. M., Tsarkova, A. S., & Yampolsky, I. V. (2016). 1001 lights: luciferins, luciferases, their mechanisms of action and applications in chemical analysis, biology and medicine. Chemical Society Reviews, 45(21), 6048-6077.
[4] Rao, J., Dragulescu-Andrasi, A., & Yao, H. (2007). Fluorescence imaging in vivo: recent advances. Current opinion in biotechnology, 18(1), 17-25.


text:
How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms

I am working with RNA-seq data without a reference genome/transcriptome and am instead using a Trinity de novo transcriptome assembly.
I analyzed both isoform and gene expression abundances using RSEM.
And got differentially expressed GENES via DESeq2.
In the next step I want to get fasta sequences of differentially expressed genes to use them in the subsequent BlastX/Blast2Go analyses.
Considering the fact that my Trinity created assembly is based on isoform id:
How can I extract the gene sequence out of it? Is that a right thing to combine all the isoforms of a single gene in order to generate the longest possible isoform assembly?
If so, how can I do it to have a union-exon (combined form of all the exons present in the isoforms of an individual gene) sequence used in my Blast analysis?
If not, how should I proceed with the data type available?
Trinity assembled file includes the isoform column as:
TRINITY_DN1780_c0_g1_i2
TRINITY_DN1780_c0_g1_i1
TRINITY_DN1780_c0_g1_i3
TRINITY_DN1780_c0_g1_i4
TRINITY_DN1780_c0_g1_i5

and the DESeq2 produces the output (gene differentials) as:
TRINITY_DN1780_c1_g1

From your question it isn't clear if you are trying to identify the isoforms that have a very high or very low differential expression value, or if you are trying to correctly find an orthlog for your top hits using BLAST.
If your question is about isoform-level quantification, identifying DE at the isoform level is challenging and DESeq2 is only designed for gene-level analyses. That is why it is only outputting genes rather than a specific isoform.

Many thanks for your quick reply. I attempt to identify all possible differentially expressed transcripts (resulted from comparing the samples of tw) by Blast or Blast2Go approaches (since I have no reference genome/transcriptome in hand). If at gene level, my first problem is that I can't get the fasta files of gene level (DESeq2 output) from Trinity fasta file. Is there any automated software/command to get gene level sequence from Trinity fasta file? 2. If not, my second problem is that which procedure should I follow to analyze my data at isoform level (preferably for those with high diffe

So, just to be clear - your final goal is to find the specific isoforms with the highest differential expression from your transcriptome?

There's no clear way to go from a trinity "gene" (this is just a cluster of isoforms with shared sequence) to a linear gene sequence that will be blastable and not cause a number of headaches (e.g., due to swapping the order of exons).
What I would suggest that you do is either Blast2GO the isoforms within the Trinity-defined "genes" or take the highest expressed isoform within a "gene" as the sequence you feed into Blast2GO for that gene. Alternatively just use isoforms for everything and avoid this whole issue.

There are a lot of questions in this post, and the main ones have been answered.  Here is an answer to the one that has not been answered yet:
How can I do it to have a union-exon (combined form of all the exons present in the isoforms of an individual gene) sequence used in my Blast analysis?
I don't think you should do this, but if you want one of the usually many plausible linearized representations of transcript isoforms, then Lace would do the trick.

It appears from the comment in another question that this question is also about transcript-level differential expression analysis.
If instead of using RSEM you use kallisto for transcript quantification, then you could use sleuth to perform transcript-level differential expression analysis.


text:
ncbi...eutils/esearch.fcgi: How to search now in the snp database?
According to the 2017(!) textbook, it worked at that time, for instance,
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=autism
At this time, this does not work and for each conceivable term it will find 0 records.
What is the right way to search in snp, valid in December 2018?

Welcome to the site Slepecky. Could you provide a reference to the textbook you mention? Maybe there is an updated version someone can look up. Have you tried looking the eutils manuals for this task?

Also, please explain what you are trying to do so we can come up with alternatives.

eUtilities for dbSNP is in the process of being retired. See this post from NCBI. You may want to check out the new API.
From the NCBI Insights blog post:

Entrez search is available only for human. We will retire all Eutils report formats (RS docsum XML and ASN.1, FASTA, and FLAT text) later this year. Only eSummary XML report will be supported.
Batch Query Service for all organisms will be retired October 2018. Users of human dbSNP data should transition to the new API.
We will place non-human FTP data in the archive directory. Please see the phase-out plan announced last year.


text:
Parsing gtf file for transcript ID and transcript name
I used this previously to get ensembl ID and its gene name from the gencode gtf annotation file
This was the script I used
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
        printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt

And the input file:
##description: evidence-based annotation of the human genome (GRCh38), version 21 (Ensembl 77)
##provider: GENCODE
##contact: gencode@sanger.ac.uk
##format: gtf
##date: 2014-09-29
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  transcript  12010   13670   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12010   12057   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 1; exon_id "ENSE00001948541.1"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12179   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 2; exon_id "ENSE00001671638.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12613   12697   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 3; exon_id "ENSE00001758273.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12975   13052   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 4; exon_id "ENSE00001799933.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";

What I'm trying to do now is replace the gene_ID with transcript_ID and the gene_name with transcript_name but I'm doing something wrong as I'm not getting the transcript_ID but the gene Ensembl ID instead. I'm not sure if I'm choosing the right column for parsing the transcript_ID and name.
So I'm expecting to get this
transcript_id "ENST00000450305.2";  gene_name "DDX11L1";transcript_name "DDX11L1-001"

Any help or suggestion would be highly appreciated

I get the Ensembl transcript ID and the transcript name (e.g., "ENST00000320153 DDI2-201"). What is the output you're getting?

im getting ENSG and ENST  ids both along with gene name, did you run the same thing?im doing something wrong perhaps

I'm running the exact same thing with /transcript_id|transcript_name/ rather than /gene_id|gene_name/.

okay ,yes im doing the same thing...but my output is everything combined gene as well as transcript ID with names if you can look into my file may be you can tell https://drive.google.com/open?id=1Frzy0R1j1Vncog6LP5E8N2TIJFBZsS3z what i need is can i get only transcript ID and name...

I don't see how your awk command would work since you're using whitespace as a field delimiter. In any case, you can use a much shorter perl command for this:
$ perl -lne 'print "@m" if @m=(/((?:transcript_id|gene_name|transcript_name)\s+\S+)/g);' Homo_sapiens.GRCh37.70.gtf
transcript_id "ENSG00000223972.5"; gene_name "DDX11L1"; transcript_name "DDX11L1";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";

The scriptlet will print the contents of the @m array if this line matches what you need. The regex magic is:

(?: ): this is a non-capturing group. We need this so that @m only holds what is inside the capturing parentheses (see below).
(  (?:transcript_id|gene_name|transcript_name)\s+\S+) : the outer parentheses will capture everything matched: one of transcript_id|gene_name|transcript_name followed by 1 or more whitespace characters (\s+) and one or more non-whitespace characters.
print @m if @m=... : only print if the match was successful.

i was in hurry and dire need almost all the work i do it in R ,i tried it in R due to the size of the gtf it just froze ,  looks like I need to learn some basics of awk along with perl that will make my life simpler...thank you very much i will use your code...

that was super fast in takes time with awk but this is the time it took
real 0m16.112s
user 0m14.653s
sys 0m0.693s


text:
Understanding the initialization of the Needleman-Wunsch algorithm
I would like to ask why when we initialize 0 as the starting point, there are 3 gray boxes that are not used as shown (highlighted in the red box) in the pic?

 Are these boxes to show that they are gaps or just to align them properly?

I would also like to ask how dynamic programming plays a part in filling up the rows and columns in the matrix. (I'm new to bioinformatics)

Welcome to the site Heeh. You can [edit] your question to add things to the question like more information. Could you clarify what do you mean by "how dynamic programming plays a part"? Do you want to know what is dynamic programming? Or perhaps why it is needed for this algorithm?

Hi I would like to know both about dynamic programming and why it is needed for this algorithm. I searched the meaning of dynamic programming but I do not understand the meaning behind recursive.

The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment
The dynamic programming is useful because you need first to calculate the score for each position, and then once you have the score, you can calculate the optimum path. So you solve your problem by first calculating what would happen if (there was a gap, there was a mismatch...) and then you find the optimum alignment.

Dynamic programming isn't strictly needed, it's just convenient to use recursion to handle the rows (depending on the implementation, you can also save memory by not storing the full matrix then).


text:
SingleCellExperiment: invalid rownames length
I did all of the installations from scratch in a virtual conda environment. Then I tried to run R script from within javascript in the following way:
let shCMD = 'Rscript --vanilla R/full_analysis.R ';
  shCMD += '--diff_1 ' + diff_str_1 + '--diff_2 ' + diff_str_2;
  shCMD += '--markers ' + filterStr + '--markers-highlighted ' + colorStr;
  shCMD += '--datasets ' + datasetStr;

  console.log('constructed command for full analysis: ');
  console.log(shCMD);
  shell.exec('source activate rmain', function(env_code, env_stdout,
env_stderr) {
    console.log('activated rmain environment')
    if (env_code === 0) {
      shell.exec(shCMD, function(code, stdout, stderr) {

As can be seen I am activating the environment by source activate rmain and I am pretty sure that the next shell command that invokes R script runs in that environment because R script needs argparse package and it is installed only in the rmain virtual environment: in rmain argparse package is successfully imported, but outside of it it gives an error for the argparse absence.
So, now the following line in the R script:
sce <- SingleCellExperiment(list(counts=as.matrix(data[,-1])))

gives me the error:
Error in `rownames<-`(`*tmp*`, value = c("0610007P14Rik", "0610009B22Rik",  :
  invalid rownames length
Calls: SingleCellExperiment ... elementMetadata -> .local -> rownames<- -> rownames<-
Execution halted

However, if I launch Rstudio from within the rmain environment, the line runs successfully... So, I am puzzled now on what is going on.

If it works interactively (in Rstudio) but not as a script it could be because R loads certain libraries for you in an interactive session. You could try adding library(methods) to the script and see if you still have the same problem


text:
Extracting all reads from bam file which match read IDs in another file
I have a long list of read IDs of interest to me in a file called read_names.txt. it is simply in the format:
m54197_180831_211346/4981510/ccs
m54197_180831_211346/6226723/ccs
m54197_180831_211346/6619607/ccs
...

etc where these are the actual read names from a fastq file. I am then trying to find those in a bam file for which the fastq has been mapped, and I can accomplish this like:
for ID in `cat read_names.txt`
do
samtools view inbam.bam | grep $ID >> read_locs.sam
done

However, this method is obscenely slow because it is rerunning samtools view for every ID iteration (several hours now for 600 read IDs), and I was hoping to do this for several read_names.txt files.
I tried sort of flipping the script a bit and running samtools view first but it only returned the first read ID present in the file and stopped:
samtools view inbam.bam | for ID in `cat read_names.txt`; do grep $ID >> read_locs.txt; done

Am I missing something or is this the only way to accomplish this task?

Just as a general rule never do for var in $(cat file). Also known as bash pitfall number 1.

thanks terdon, I will keep this in mind and find the better ways to handle this in the future. I use this method quite often so that is great for me to know!

Yeah, most of us use it very often until it bites us! :)

It is still slow but grep has a -f option to take in a file
samtools view inbam.bam | grep -f read_names.txt > read_locs.txt

This took 30 seconds as opposed to several hours, so I'd say the overhead was definitely in samtools view. I feel silly because of how simple this was. Thank you!

just as a warning this can take much longer if the read name file is big.

I appreciate the heads up! Comparatively, I do not think it will ever take longer than the method I was using so it is still a great answer. (This is serving as a caveat to your warning).

samtools can do this natively too using -N, --qname-file CLI option:
samtools view -N read_names.txt in.bam > read_locs.sam

❯ samtools --version
samtools 1.16.1
Using htslib 1.16