text:
Problem with classification model of genomic data: every machine learing model predicts wrongly almost always the same subset of dataset
First of all, I'd like to apologize for any spelling or grammar mistakes.
I'm having a problem using R for a classification problem. My dataset contains ~300.000 genomic data, and the features are DNA-related features (number of dinucleotides, number of trinucleotides, the CG Content, and some more). In conclusion, I have a dataset of 300.000 rows and 84 columns (columns= features). The 84th feature is basically the classification variable (there are two classes: class 1 and class 2).
I split the dataset in train set and test set and fit the train set to create a model using different methods, like Random Forest, XGBoost, Adaboost, etc. and check the accuracy using the test set.
In every model, the accuracy of the prediction on the test set is always around 86-87%. What I noticed is that, in the test set, for EVERY MODEL, almost the same number of rows (data) is predicted right (the model predict its true class), and always almost the same data of rows (data) is predicted wrong (the model predict wrong class. For example its true class is 1 but the model predicts 2).
If I take the set of data that are predicted right for every model I've created (250.000 data from the 300.000), and build a new model out of them (splitting them again in a new train and test set), I have an accuracy of 98%.
That means that the other 50.000 data (from the 300.000), add an error at every model I'm trying to create.
But I can't just ignore them. I thought that it would be a way to create two differents models, but doesn't know how to figure this out.
Then I thought that there could be a way to find which variables from the two subsets have strong differences, but even if I figure this out, I don't know what's the way for example to scale the original dataset to prevent those differences from keeping the accuracy in 87% in every model (using the original dataset of the 300.000 rows).
In conclusion, I have a dataset that pretends like "two datasets". It consists of two subsets: the first one is always predicted right from any model (machine learning algorithm) I'm trying, and the second one is always predicted wrong. How could I build a strong model taking into account the differences?
I also took a look at the box plots of the features. The most important feature for every model is called "CG" (number of dinucleotide of cytosines and guanines in DNA). In the following image you can see the boxplot of CG feature separating the original dataset in two datasets (T = for the dataset that is predicted right in every model) and (F = for the dataset that is predicted wrong in every model). 1 stands for data of class 1 and 2 stands for data of class 2. So, T1 means data from the right predicted dataset of class 1, T2 means data from the right predicted dataset of classes 2, T means data from the right predicted dataset for both class 1 and 2, etc.

As I'm a new beginner in R and in machine learning, I hope you do understand my problem. Thanks sincerely.

It needs some thought. Whats your training/test split, 80:20? Any parameter splits?

The persistent 87-88% accuracy between models is a bit weird, but its an interesting problem

I always used a train set with same amount of class 1 and class 2. for example, a train set of 50.000 data (25.000 of class 1 and 25.000 of class 2).

Okay thats the classificaiton (I don't know what "Class I" and "Class 2" represent, but I get the idea in ML, whatever "Class" is, is the thing you are looking at. The question was the split with testing, or parameterisation. The reliability of the accuracy score across very different models suggests something is not right, but I honestly don't know what.

Are your classes balanced, meaning that do you have roughly the same number of the two classes? Looks like your models are learning a particular set of samples very well at the expense of learning others. Cross-validation could help here.

Not sure it is related but do you shuffle your datasets before training ? I had weird patterns of valid accuracy from epoch to epoch when training a neural network once, it was because I forgot to shuffle my datasets.

You could be seeing the effect of SNPs on haplotype-driven conditions (i.e. where a combination of multiple variants predicts the class, rather than a single variant).
Do you have any method to determine the age of particular variants (e.g. via an LD plot to find inconsistencies between variants and the surrounding block structure), so that you can see if that correlates with the likelihood of error?

No, unforetunately no.


text:
Finding discriminating positions for a sequence alignment column
I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.
So, here is my code for the Monte Carlo method:
def MonteCarlo(S, max_iter, threshold):
    """
    Computes the discriminating position for a sequence column

    Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
    Output: 2 classes
    """
    n = len(S)
    S_copy = S.copy()
    num_class1 = int(0.30*n) #number in class1
    print(num_class1)
    for _ in range(num_class1): #intialisation
        rand_position = random.randint(0,len(S)-1)
        if S_copy[rand_position] == 'X':
            rand_position = random.randint(0,len(S)) #reselect if used position
        S_copy[rand_position] = 'X'

    #Monte Carlo with Insertions and Deletions

    s_entropy = entropy(S)
    n_iter = 0 #counter for iterations
    class_1 = dict() #class 1 determined by X's
    class_2 = dict() #class 2 is rest of S
    for i in range(len(S)):
        if S_copy[i] == 'X':
            class_1[i] =S[i]
        else:
            class_2[i] =S[i]

    res = score(S, class_1, class_2)

    while n_iter < max_iter:
        r = random.random()

        if r < 0.25: #25% chance to do an insertion
          (c1,c2,pos) = Insertion(S, class_1, class_2)

        elif r < 0.5: #25% chance to do a deletion
          (c1,c2,pos) = Deletion(S, class_1, class_2)

        else: # 50% chance to do a swap
          (c1,c2,pos) = Swap(S, class_1, class_2)

        res2 = score(S, c1, c2)

        #Boltzmann Distribution coefficient
        delta_S = abs(res2 - res)
        lambda_, N = 0.2, 19
        size = 6000
        ind1 = random.randint(0,size - 1)
        rv = boltzmann.rvs(lambda_, N, size=size)/N
        beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
        if  (res2 > res) or (beta*delta_S > threshold): #if better score i.e. higher entropy gain
            res = res2
            class_1 = c1
            class_2 = c2

        n_iter += 1

    return class_1, class_2, pos

Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)
Let me explain what my code is trying to do: Given a sequence alignment (as a list) say ['M','M','M',...,'L','L','L'] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.
While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results.
Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?
Thanks!

It's not very clear what you mean by "discriminating position". Can you explain more or link to a source that describes in a little more detail what you are talking about?

Hi! By discriminating position, I mean I'm trying to find a position within the sequence alignment that splits the columns into two classes such that the two have the lowest entropy possible. So, in my sample case with a clear division between Ms and Ls, I should've gotten the the position at which the Ls begin. Does that make sense?
Here is a link: https://pubs.rsc.org/en/content/articlelanding/2019/sc/c8sc03669a#!divAbstract

Hi, the paper that you've linked to is talking about a discrimination between two structural states of a protein by a kinase; "discriminating position" does not appear to occur in its full text. It's not clear how this generalizes to alignments. From your explanation though I think I have a better sense of what you mean: for seq alignment $S$ with positions i = 1,2,3...n find the position $S_i$ such that a measure of entropy across aligned positions $S_1 : S_{i-1}$ and $S_{i+1} : S_n$ is minimized for both slices across each column of the aligned slice. (words are easier than code for me).

Yes! no problem (: I also misunderstood what I was trying to solve actually, I have managed to fix the problem! But for future reference: the code is supposed to split the sequence column into two classes such that each distinctly has the highest entropy possible (or lowest for that matter)

Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.

Right I see, so I need to preprocess the sequence alignments first? So that I match the gaps properly?

Yes, to verify this step for logical errors

source: Stack Exchange

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

sort: votes


text:
How to analyze co-occurrence of multiple SNPs?
I am interested in 20 different SNPs that all are either As or Gs, and they all occur on the same chromosome. How can I assess the co-occurrence of these SNPs? In other words, I want to know, if SNP1 is "G", does that make SNP2 to be more likely to be a "G" as well?
I know this is broad and I'm not looking for a script or anything so specific, just a point in the right direction. I'm working with a non-model organism, so any organism-specific application won't be useful without being able to learn the methodology.
Googling thus far has only brought up GWAS and similar methods of correlating a genotype and a phenotype, but not within a genotype.

Okay, just found the term "linkage disequilibirum" which should be helpful...

I agree with @CephBirk that you are just looking for linkage disequilibrium. To calculate this in a simple way, it is just the Pearson correlation coefficient squared.

I found the LDheatmap package in R very helpful at calculating linkage disequilibrium. You can make plots such as this:

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


text:
How do I pull singe cell RNA sequencing data from GEO database?
I am new to R and computational biology. I am trying to look through a published data set to check for gene expression for my own project. I am having trouble finding materials to teach me how to analyze single cell RNA sequencing data on R from the GEO NCBI database.
This is the data that I am trying to access https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138651
I would appreciate any directions or suggestions for tutorials.
Thanks!

Install package
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

Use function getGEO to download the files.
expressionSet <- getGEO("GSE138651")
expressionSet <- expressionSet[[1]]

You can then look at the info of the samples:
pData(expressionSet)

And the assay data itself
exprs(expressionSet)

This is useful because of the very specific instructions but I think it would be even more helpful to also include a brief introductory explanation and links to BioConductor and the GEOquery package to provide context for newcomers.

I am sorry but I have to flag this for removal since it's wrong. getGEO only returns expression data for arrays, not for NGS. Output of exprs(expressionSet) is empty and this is expected.

The simplest would be using a count matrix (at the end of the link you have shared, section "Supplementary file"). For example GSE138651_Vagal_WholeNodoseSeq_raw.mtx.gz corresponds to "raw" count matrices obtained after mapping and counting the reads. This file, along with the GSE138651_barcodes.tsv.gz and GSE138651_genes.tsv.gz can be used to "read" the data into R; to that end there are many options in the form of packages, Seurat is one of the best documented and hence is beginner-friendly.

Thanks Haci. I am learning Seurat as well. So many things!!
I am actually the most interested in the file GSE138651_NGTall_Reads_cap.csv.gz. Do you know if this file type can be "read" into R?

Yes, it is a csv file and once again there are many options to read such files into R. read.csv() is one option. Personally, I use fread() of data.table to read in tabular files.

I will try that then. Thank you again!

You might want to check out the BioConductor GEOquery package.
Install using
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

More information on how to use this package and read in data from GEO can be found in the vignette

I am sorry but I have to flag this for removal since it's wrong. getGEO only returns expression data for arrays, not for NGS. This package does not apply here.


text:
Is it possible to get such a plot
I wish to reproduce a very informative plot in google

Given a group of samples and a list of genes for which we have data how each genes has mutated and whether the alleles are homozygous of heterozygous. These are expressed as boolean values
Genes enriched in 2 pathways like cell cycle and p53 as shown in this figure
Map of functional alterations for a group of patients. Genes (rows) encoding components p53–DNA repair pathway;  are affected by selected functional events (percent of samples altered and types of alteration are represented by colored squares) across group of samples in column. Alterations of the pathway are observed stacked green bar plots at bottom
I have the boolean matrix of mutated genes per sample would be grateful for assistence.
An example of my copy number data
> dataset_gistic$annotations
     type                event
G2   "Homozygous Deletion"   "ACVR2A"
G13  "Homozygous Deletion"   "HIST1H3B"
G73  "Heterozygous Deletion" "CCNE1"
G74  "Heterozygous Deletion" "TSHZ3"
G78  "Low-level Amplification"    "CTNNB1"
G110 "Low-level Amplification"    "CCNE1"
G111 "Low-level Amplification"    "TSHZ3"
G119 "High-level Amplification"   "PIK3CA"
G121 "High-level Amplification"   "NIPBL"
G147 "High-level Amplification"   "CCNE1"
G148 "High-level Amplification"   "TSHZ3"
>

> dput(dataset_gistic$annotations)
structure(c("Homozygous Loss", "Homozygous Loss", "Homozygous Loss",
"Homozygous Loss", "Homozygous Loss", "Homozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "ACVR2A", "HIST1H3B", "CDKN2A", "PTEN", "NAV3",
"LIN7A", "ARID1A", "ACVR2A", "SCN3A", "CTNNB1", "PBRM1", "EPHA3",
"POLQ", "PIK3CA", "SLIT2", "NIPBL", "MSH3", "APC", "HIST1H3B",
"ARID1B", "CDKN2A", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "NAV3", "LIN7A", "SIN3A", "TP53", "ERBB2", "GPATCH8",
"RNF43", "GATA6", "SMAD4", "CCDC102B", "STK11", "SMARCA4", "CCNE1",
"TSHZ3", "CTNNB1", "PIK3CA", "SLIT2", "NIPBL", "HIST1H3B", "ARID1B",
"EGFR", "KCNQ3", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "SIN3A", "ERBB2", "GPATCH8", "RNF43", "GATA6", "CCDC102B",
"STK11", "SMARCA4", "CCNE1", "TSHZ3", "PIK3CA", "NIPBL", "EGFR",
"KCNQ3", "CCND1", "KRAS", "ERBB2", "GPATCH8", "RNF43", "GATA6",
"CCNE1", "TSHZ3"), .Dim = c(77L, 2L), .Dimnames = list(c("G2",
"G13", "G17", "G19", "G24", "G25", "G38", "G39", "G40", "G41",
"G42", "G43", "G44", "G45", "G46", "G47", "G48", "G49", "G50",
"G51", "G54", "G55", "G56", "G57", "G58", "G59", "G60", "G61",
"G62", "G63", "G64", "G65", "G66", "G67", "G68", "G69", "G70",
"G71", "G72", "G73", "G74", "G78", "G82", "G83", "G84", "G87",
"G88", "G89", "G90", "G92", "G93", "G94", "G95", "G96", "G97",
"G100", "G102", "G103", "G104", "G105", "G107", "G108", "G109",
"G110", "G111", "G119", "G121", "G126", "G127", "G131", "G132",
"G139", "G140", "G141", "G142", "G147", "G148"), c("type", "event"
)))
>

And this is sample of my mutational data the number of mutations for each gene in each sample
    om$numericMatrix
       LP6005500-DNA_D03 LP6008336-DNA_G01 LP6008460-DNA_D01 LP6005334-DNA_H01
TP53                   2                 2                 2                 2
CDKN2A                 3                 1                 2                 0
       LP6007600 LP6008334-DNA_A03 LP6008334-DNA_A04 LP6008334-DNA_B02
TP53           2                 2                 2                 2
CDKN2A         0                 0                 0                 0
       LP6008334-DNA_C02 LP6008334-DNA_D02 LP6008336-DNA_H01 LP6008337-DNA_A07
TP53                   2                 1                 2                 2
CDKN2A                 0                 0                 0                 0
       LP6008337-DNA_H06 LP6008460-DNA_A04
TP53                   2                 2
CDKN2A                 0                 0
>

> dput(om$numericMatrix)
structure(c(2, 3, 2, 1, 2, 2, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2,
0, 1, 0, 2, 0, 2, 0, 2, 0, 2, 0), .Dim = c(2L, 14L), .Dimnames = list(
    c("TP53", "CDKN2A"), c("LP6005500-DNA_D03", "LP6008336-DNA_G01",
    "LP6008460-DNA_D01", "LP6005334-DNA_H01", "LP6007600", "LP6008334-DNA_A03",
    "LP6008334-DNA_A04", "LP6008334-DNA_B02", "LP6008334-DNA_C02",
    "LP6008334-DNA_D02", "LP6008336-DNA_H01", "LP6008337-DNA_A07",
    "LP6008337-DNA_H06", "LP6008460-DNA_A04")))
>

Why unvoting? Do you think only by a simple googling one can get such a plot?

Can you update the question by giving an example of your data looks like instead of describing them by words?

@Exhausted, please provide the community with a reproducible example: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example As was suggested in another question of yours, dput() will help us "reproduce" your input data. If you invest a couple of minutes going through the link and subsequently increase the "quality" and "completeness" of your questions, you will get answers much faster you do now.

While it was easy to get your data after your edit, the data was not quite useful, the second data frame/matrix has 2 columns only!
I had to create a fake data set, that looks like this:
my_df <- data.frame(
  gene_class = rep(letters[1:4],3),
  some_value = rep(1:3*10,4)
)
rownames(my_df) <- paste("Gene-",rep(LETTERS[1:4],each = 3),1:4,sep="")

> my_df
        gene_class some_value
Gene-A1          a         10
Gene-A2          b         20
Gene-A3          c         30
Gene-B4          d         10
Gene-B1          a         20
Gene-B2          b         30
Gene-C3          c         10
Gene-C4          d         20
Gene-C1          a         30
Gene-D2          b         10
Gene-D3          c         20
Gene-D4          d         30

Using the ComplexHeatmap package, I got something similar to the figure you have posted.
# data frame for annotation
my_annot <- my_df[,"gene_class",drop=FALSE]

# the Heatmap() function requires a matrix
my_matrix <- as.matrix(my_df[,"some_value",drop=FALSE])

dev.off() # clear the plotting area in case not empty
ha = rowAnnotation(gene_class = my_annot$gene_class)

ht <-  Heatmap(my_matrix,
               left_annotation = ha, # gene_class
               row_names_side = "left", # position of row (gene) names
               heatmap_width = unit(0.3, "cm")*nrow(my_matrix), # cell width
               cluster_rows = FALSE, # no clustering
               cluster_columns = FALSE, # no clustering
               cell_fun = function(j, i, x, y, w, h, col) { # add text ("some value column") to each grid
                 grid.text(my_matrix[i, j], x, y)
               })
draw(ht, heatmap_legend_side = "bottom")

The ComplexHeatmap package offers enormous flexibility for customization. I strongly recommend devoting some time to it. The way I use it, I try to put all of my data/information in a single data frame/table, which I subsequently subset to create data frames and matrices for annotations and the actual heatmap(s) respectively.

Thank you @haci . You are amazing and helpful as always, definitely your solution works but I don't know how to format my own data to work with your code. What you have produced is what I need as long as I am able to adapt my data our your dummy input


text:
Output when compute GC content is different between using user input and using string
I am learning bioinformatics and i was solving a problem on rosalind.info.
Input is a string consist of multiple DNA strand names and the DNA sequences, and Output is the name of the sequence that has the highest GC percentage and the GC percentage.
this is the input:
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

this is the correct output:
Rosalind_0808
60.919540

This is my code (sorry for this filthy mess):
DNA = input()
wholeRosalind_DNA = DNA.splitlines()
Rosalind = []
Rosalind_DNA = []
for i in wholeRosalind_DNA:
    if "Rosalind" in i:
        Rosalind.append(i)
Countstring = []
maxGCcontent = []
GCcontent = []
def eachRosalind_DNA(DNA):
    for i in range(0,len(Rosalind)):
        if Rosalind[i] in wholeRosalind_DNA:
            a = wholeRosalind_DNA.index(Rosalind[i])
            if a == wholeRosalind_DNA.index(Rosalind[-1]):
                b = len(wholeRosalind_DNA)+1
            else:
                b = wholeRosalind_DNA.index(Rosalind[i+1])
            Rosalind_DNA.append(wholeRosalind_DNA[a:b])
    for i in Rosalind_DNA:
        Countstring.append(''.join(i))
    for i in Countstring:
        bases = len(i)-14
        CG = i.count("C") + i.count("G")
        GCcontent.append([i[1:14],(100*(CG/bases))])
    for i in GCcontent:
        maxGCcontent.append(i[1])
    result = max(maxGCcontent)
    for i in GCcontent:
        if i[1] == result:
            return (i[0]),(round(i[1],6))
print(eachRosalind_DNA(DNA)[0])
print(eachRosalind_DNA(DNA)[1])

but i get this instead:
  File "GC_Content.py", line 26, in eachRosalind_DNA
    GCcontent.append([i[1:14],(100*(CG/bases))])
ZeroDivisionError: division by zero

i will get the desired output if i just paste the input into the variabe "DNA" without using user input. Please help me fix the code and distinguish those two ways of getting input

Piece of advice: use meaningful names for your variables (loop variables do not need to be named i). This will help you and others spotting potential bugs.

Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:
from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

Output:
('Rosalind_0808', 60.91954022988506)

The error is that bases = 0. Briefly looking at the code carriage returns in your original code will trigger this, as will confusing a faster header with sequence. The latter is because the fasta headers are 14 characters long and subtracting 14 from bases which is of course zero. The former will likely result in 'nan' because it would be 14 subtracted from ''.
I've a simpler solution, download Biopython (using conda or last resort pip), you might want to use SeqIO (Biopython), which would be cooler because it allows you to call Rosalind_0808 directly. Then:
from Bio.SeqUtils import GC
seq = 'CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT'
print (GC(seq))

Done. If you want to splice off the last 14 characters its
GC(seq[:-14])

I think the hard thing in this problem is not the calculation of GC content, it is the processing of DNA strings because they are separated on different lines. I would suggest reading the file line by line to generate complete DNA strings for each Rosalind entry. (I am not able to spot what is causing the error in the code).
seqDic = {}
with open("YouInFile", "r") as inFile:
    for ln in inFile:
        ln = ln.rstrip()
        if ln.startswith(">"):
            ln.name = ln[1:]
            seq = inFile.readline().rstrip() + inFile.readline().rstrip()
            seqDic[ln.name] = seq
        else
            continue

Basically, whenever you encounter a ">" in your file, you write down the Rosalind name, and then you read the next two lines and combine them together. You put the name and the sequence in a dictionary, and iterate over the dictionary to calculate GC content.
Note that I did not debug my code. Let me know if it is not working.


text:
BigWig circular plot representation
I would like to ask if anyone knows or has experience in representing BigWig signals in a circular plot? I mean like IGV but in circular plot not linear. Or any library that can be used in plotting BigWig signals especially in Javascript?
I appreciate any suggestion,
Thank you.

circosJS should probably work well, based on your requirements: https://github.com/nicgirault/circosJS
I posted a short Observable notebook that shows what calls could be used to bring data from remote bigWig files into a circosJS plot:
https://observablehq.com/@alexpreynolds/rendering-bigwig-file-data-in-a-circular-representation
It is less intended as an explicit how-to in Javascript, than as a starting point for the end user with some of the function calls and data scale issues involved. Hopefully it can be useful to the reader for exploratory purposes.

Does it take bigwig file without any conversion to tabulated format?

Doesn't look like it, but you could probably use https://github.com/GMOD/bbi-js as a file parser, get the data you need asynchronously, and feed it to circosJS to render.

Thank you very much on your suggestion, I will give it a try and see if it works.


text:
Database of RNA-Seq and SNP genotyping data
I need to obtain publicly available RNA-seq + SNP genotyping data for a cohort - ideally one where the RNA-seq data has been sequenced from brain tissue.
I have found CMC and GTEx to be sufficient for my needs, however the procedure to obtaining data from both these sources are extremely long and they request us to share our data with them as well; which we can not do because it is not ours to share.
I was wondering if anyone knows of another resource that can be used for obtaining such data?
Thanks

Have a look through EGA. It's not public, but you're unlikely to find anything that is.

Human RNA-seq data, like genetic data is very sensitive data as it can be used to identify people (and their relatives), so people will likely not distribute it without entering a formal agreement that ensures that people can be held accountable if the data is misused, stolen or made public. Depending on the consent form signed by the patients (or their relatives), they may not be allowed to share it at all. One possibility would be to start a collaboration with a group that has access to such data and ask them to run your analyses on their data. You can then have a shared publication and it's a win for everybody.

Could you elaborate on your answer, because at present it is not sufficiently informative? "will likely not distribute it lightly", what does that mean precisely by means of example?

In most western countries making such data publicly available would be illegal, or at the very least be deemed so ethically problematic that data would be sent to EGA or a similar semi-closed platform.

Thank you for your comments. I will take a look at EGA. @citronrose very good suggestion, I thought about it myself. Unfortunately with the current climate, and with everybody working from home, it's not as easy to contact people as it once was, I will give it a go though.


text:
How I reformat this data
I have this data all is one line even gene from first sample joined the second one while I want this data column wise each field separated like I have provided
#Uploaded variation   #Uploaded_variation 1 883477 G T 1_883477_G/TPRDM16 1 3125744 A G 1_3125744_A/GPRDM16

But I want
         V1     V2  V3      V4  V5  V6
1_883477_G/T    1   883477  G   T   PRDM16
1_3125744_A/G   1   3125744 A   G   PRDM16

Hi @Exhausted, this a similar style of question to your Google heatmap question, in that you posted the input and output you wanted and then requested all the intervening code. Could post the code you have tried so far in line with the sites requirements? Please note, I didn't downvote you, but whoever did downvoted you because of your question format. In summary, Bioinfo SE is not a code writing service.

Are you sure it's one line? If you're writing a shell script, make sure you quote the variable: echo "$data"  and not echo $data

The reason I imagine you constantly get downvoted is your attitude. You have the habit of outsourcing your work on a daily basis to the bioinformatics online communities, e.g. Biostars and SE. Day after day you post the same types of questions: How to do this plot/figure, how to rearrange this dataset. This is not a code factors here. Most of the time you do not invest effort into showing what you tried. Often the posts do not even include proper data examples. It is like "here is my problem, give my solution". The idea of a community is to help each other learning things, not to do your work.

Yes you right :( I am really not happy with this even I am suffering from my attitude than the community do

I agree with the other commentators, but wanted to add that only the secondary purpose of Stack Exchange sites is to help the question asker. The primary purpose is actually to create a useful repository of questions and answers for the whole community. Unfortunately, your questions are not useful for others and titles like "How I format this data" are way too vague. Don't worry though, you can learn from this and ask better questions in the future :)

I am really curious what software or piece of code gave you this one-line data. I would think twice before using a tool giving such output. Anyway, I think regular expressions is the way to go in this case (this solution uses R):
library(stringr)

text <- "#Uploaded variation   #Uploaded_variation 1 883477 G T 1_883477_G/TPRDM16 1 3125744 A G 1_3125744_A/GPRDM16"

my_pattern <- "[0-9]+\\s[0-9]+\\s[ATGC]+\\s[ATGC]+\\s[0-9]+_[0-9]+_[ATGC]+\\/[0-9A-Z]+"

info_as_list <- str_extract_all(text,
                                pattern = my_pattern) %>% unlist()

info_as_list <- sub("/",
                    " ",
                    info_as_list)

info_as_list <- str_split(info_as_list, " ")

as.data.frame(Reduce(rbind, info_as_list))

This is the output, I guess the rowname "init" comes from Reduce():
     V1      V2 V3 V4          V5      V6
init  1  883477  G  T  1_883477_G TPRDM16
      1 3125744  A  G 1_3125744_A GPRDM16

Here is what the regular expression does, you might need to alter this based on your data:
[0-9]+\\s       :   start position
[0-9]+\\s       :   end position
[ATGC]+\\s      :   ref nucleotide(s)
[ATGC]+\\s      :   alternate nucleotide(s)
[0-9]+          :   start
_               :   underscore
[0-9]+          :   end
_               :   underscore
[ATGC]+         :   ref
\\/             :   slash
[0-9A-Z]+       :   gene symbol

+ means one or more of the preceding and R uses two backslashes, \\, as escape characters.

Sorry @haci everything is fine except for your are adding alternative allele to the start of gene like TPRDM16 while in column V5 I need 1_883477_G/T and in column V6 I need gene name which is PRDM16.

You can alter column 5 by pasting it with column 4. Use paste() with sep="/".

I don't want to paste column V5 to V6, they are independent columns. In column V6 I need only PRDM16 because this is a gene symbol and in column V5 I need a combination of columns V1 to V4 as 1_883477_G/T. Your kindly code pastes a part of column V5 to V6 which is undesired. If I paste column V5 to V4 I also have to remove the first character of column V6 which is too confusing for me

@Exhausted, this is not a code writing service as @Michael stated above. We cannot give you completely working code but can only point you to the right direction. If you would have googled "R remove the first character of a string", you would find out substring() is the way to go to "clean" the extra character in your "gene name" column: df$V6 <- substring(df$V6,2)


text:
How can I collapse probes to single value per gene using GEOquery?
I have extracted microarray gene expression data from GEO using the GEOquery package in R. There are multiple probes per gene. Is there a way to select the best probe? If so, what would be the criteria and how can you do this? Maybe use Jetset probes (=best quality probes) for this?
Otherwise, would you use the average expression value of all probes for each gene? If so, how can you do this in R, please?
EDIT
Edit after comment below:
I like the suggestion proposed here: https://support.bioconductor.org/p/70133/, using this code:
 Probesets=paste("a",1:200,sep="") # "fake probesets"
 Genes=sample(letters,200,replace=T)
 Value=rnorm(200)
 X=data.frame(Probesets,Genes,Value) # I edited this line, as it seemed to contain an error
 X=X[order(X$Value,decreasing=T),]
 Y=X[which(!duplicated(X$Genes)),]

But, I do not understand it. How is the value for each gene in Y calculated? It was commented that this is the MAD, but calculating the Median Absolute Deviation (MAD) by hand results in a different value. Could anyone help, please? Thank you.

The post: https://bioinformatics.stackexchange.com/questions/3082/duplicate-gene-symbol-handling-in-geo-gene-expression-data/3083#3083 helped a little bit, but did not answer all my questions.

https://support.bioconductor.org/p/92128/

Thank you. Lots of options! I still have a question. See edit above. Thank you.

The trick is duplicated() being used with ! and only after sorting X in decreasing order. Basically, you sort your "values" from higher to lower and duplicated() gives you FALSE if that particular "value" was never observed in the "array of interest". When you negate this with !, you only keep the first occurence of a particular value and this first occurence would be the highest as it was sorted already in the previous step.

Thank you for your answer. That is very clear now. Now I am only left with figuring out how to calculate the MAD in R, but I will search for that.


text:
Accessing .bam/.cram files from AWS S3?
What's the best/easiest way for accessing .bam/.cram files from S3?
I have .bam/.bai .cram/.crai and .bed + .gff files that I want to make available easily to others so they can browse in IGV and other solutions from their computers. There seem to be different solutions after cursory Googling, what would you recommend?
Here is what I've found so far:
IGV:

https://twitter.com/igvteam/status/1187766079457644544
https://umccr.org/blog/igv-amazon-backend-setup/

Samtools: samtools directly supports AWS

http://www.htslib.org/doc/htslib-s3-plugin.html

The full URL format is:

    s3[+SCHEME]://[ID[:SECRET[:TOKEN]]@]BUCKET/PATH

samtools directly supports AWS these days: http://www.htslib.org/doc/htslib-s3-plugin.html

You could set up a public read-only bucket in an S3 account and put your files there. URLs would look like:
https://${username}.s3.amazonaws.com/public/${filename}

For samtools, the request for a BAM file's index should be built off that URL. For IGV, I think you can specify separate, custom URLs for the data and the index resources.
Public buckets are available to all, without need for credentials. So that might be an issue to think about, with respect to billing for web requests and data transfer.
I'm still not 100% certain that samtools/tabix support https requests reliably (I have a prior question on here about that), but that might be one approach.

Doo you have your heart set on AWS, or are you open to using Google Cloud Platform (GCP)?
If you're open to using GCP, there are some advantages. For example GATK can read and write directly to and from Google storage buckets, and IGV can use your login credentials to access data in Google storage as described here:
https://googlegenomics.readthedocs.io/en/latest/use_cases/browse_genomic_data/igv.html
Note that if you make the data public, by default you're going to get billed if people download it. Not sure how that's handled on AWS, but on GCP you can set the bucket (or just the file) to be "requester-pays", ie the person who downloads is the person who pays for it:
https://cloud.google.com/storage/docs/requester-pays
In general there are several ways to interact with Google buckets.

web-based file browser (example is a public bucket I manage): https://console.cloud.google.com/storage/browser/genomics-in-the-cloud
command-line access through the cloud sdk toolkit gsutil: https://cloud.google.com/storage/docs/gsutil

Disclosure: I work at the Broad Institute, which uses a lot of GCP infrastructure; I do not get any compensation from Google Cloud, but I recently wrote a book about Genomics in the Cloud that uses GCP.

I have recently open sourced a Java NIO SPI for S3 that allows direct reads from S3 using Java applications that use the NIO model. This allows applications like GATK and some parts of Picard to understand s3:// URIs.
It is used as a drop in so as long as the jar is on your classpath or in your JRE ext folder then Java will delegate nio operations to this lib when it "sees" an s3:// URI. It will use the credentials and configuration of your current profile so no additional config is require and no recompilation of an application is needed.
The github repo is https://github.com/awslabs/aws-java-nio-spi-for-s3
Disclosure: I am the author of the library, employed by the AWS Health AI Genomics Team, and have financial interests in Amazon


text:
SNPEFF genome build with specific codon usage
Failing to build my own Arabidopsis thaliana reference with Mt using 'codon.Mitochondrial' and Pt using 'codon.Bacterial_and_Plant_Plastid' but the same commands works well with both lines with codon tables removed below.
The snpeff commands
java -Xms24g -Xmx24g -jar /opt/biotools/snpEff/snpEff.jar build -gff3 -v Arabidopsis_thaliana.TAIR10.46

returns:
java.lang.RuntimeException: Error parsing property 'Arabidopsis_thaliana.TAIR10.46.Pt.codonTable'. No such codon table 'codon.Bacterial_and_Plant_Plastid'
        at org.snpeff.snpEffect.Config.createCodonTables(Config.java:169)
        at org.snpeff.snpEffect.Config.readConfig(Config.java:650)
        at org.snpeff.snpEffect.Config.init(Config.java:480)
        at org.snpeff.snpEffect.Config.<init>(Config.java:117)
        at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
        at org.snpeff.SnpEff.run(SnpEff.java:1183)
        at org.snpeff.SnpEff.main(SnpEff.java:162)

(same error for the Mt when removing Pt)
my fasta has the two chromosomes
>1 dna:chromosome chromosome:TAIR10:1:1:30427671:1 REF
>2 dna:chromosome chromosome:TAIR10:2:1:19698289:1 REF
>3 dna:chromosome chromosome:TAIR10:3:1:23459830:1 REF
>4 dna:chromosome chromosome:TAIR10:4:1:18585056:1 REF
>5 dna:chromosome chromosome:TAIR10:5:1:26975502:1 REF
>Mt dna:chromosome chromosome:TAIR10:Mt:1:366924:1 REF
>Pt dna:chromosome chromosome:TAIR10:Pt:1:154478:1 REF

my data folder is OK
data/Arabidopsis_thaliana.TAIR10.46/
drwxr-xr-x   3 u0002316 domain users 4.0K Apr 15 13:33 .
drwxr-xr-x 244 u0002316 domain users  16K Apr 10 09:39 ..
-rw-r--r--   1 u0002316 domain users 106M Apr 10 09:40 genes.gff
-rw-r--r--   1 u0002316 domain users 117M Apr 10 09:40 sequences.fa

my snpeff.config contains:
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : codon.Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : codon.Bacterial_and_Plant_Plastid

the two codon tables are present by default in the snpeff.config
codon.Mitochondrial                         : TTT/F, TTC/F, TTA/L+, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/W, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

codon.Bacterial_and_Plant_Plastid           : TTT/F, TTC/F, TTA/L, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/*, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

Anybody ever succeeded to build a custom genome with specific codon tables?
Thanks in advance

codon table issue
I solved the issue with codon tables, the correct way to specify a table is the following:
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : Bacterial_and_Plant_Plastid

Note that one should remove the leading prefix 'codon.' from the name of the tables chosen from the config file
Concerning accessory files:
data to reproduce this issue:
baseurl="ftp://ftp.ensemblgenomes.org/pub/plants/release-46"
wget ${baseurl}/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ${baseurl}/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gff3.gz
wget ${baseurl}/variation/vcf/arabidopsis_thaliana/arabidopsis_thaliana_incl_consequences.vcf.gz

I suspect that the format used by Ensembl is not 100% OK for SnpEff; I therefore simplified the fasta headers to keep only the first argument (name) and remove the remaining of the line
eg.
>AT3G05780.1 pep chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 transcript:AT3G05780.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]

becomes
>AT3G05780.1

Despite this modification, I still have an issue with
java -jar /opt/biotools/snpEff/SnpEff.jar build -gff3 -v Arabidopsis_thaliana_TAIR10_46 2>&1 | tee -a log.txt

a LOT of warning about not being able to add UTR while it seems that UTR's are well present in the GFF3 file

Reading GFF3 data file  :
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  WARNING: Cannot find transcript 'transcript:AT1G01010.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Cannot find transcript 'transcript:AT1G01020.3'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Cannot find transcript 'transcript:AT1G01020.5'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Cannot find transcript 'transcript:AT1G01020.4'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Cannot find transcript 'transcript:AT1G01020.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Cannot find transcript 'transcript:AT1G01030.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line 155  '1  araport11   three_prime_UTR 11649   11863   .   -   .   Parent

nothing good coming out of the cds.fa input although I treated it with the same name fix as the pep file

00:00:19  CDS check (FASTA file):
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'
00:00:19  Checking database using CDS sequences 00:00:19  Reading CDSs
  from file
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'...
  00:00:20  done (75976 CDSs). 00:00:20 Comparing CDS...    Labels:         '+' :
  OK        '.' : Missing       '*' : Error
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
many more lines with dots
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
    ...............................     CDS
  check:    Arabidopsis_thaliana_TAIR10_46  OK: 0   Warnings: 0 Not found:
  104147    Errors: 0   Error percentage: NaN% Fatal error: No CDS checked.
  This is might be caused by differences in FASTA file transcript IDs
  respect to database's transcript's IDs. Transcript IDs from database
  (sample):     'transcript:AT3G12250.7'    'transcript:AT5G18570.1'
    'AT3G27810.1'   'AT3G27810.2'   'transcript:AT2G16810.1'
    'transcript:AT3G01170.1'    'AT4G37030.1'   'transcript:AT1G79940.4'
    'transcript:AT2G44850.1'    'AT1G49540.1'   'AT1G49540.2'
    'transcript:AT4G35290.5'    'transcript:AT4G15840.2'
    'transcript:AT4G27500.1'    'AT4G37022.1'   'transcript:AT2G40008.1'
    'AT3G27820.1'   'transcript:AT4G24265.1'    'transcript:AT2G20618.1'
    'transcript:AT4G18375.2'    'AT4G37020.2'   'AT4G37020.1' Transcript IDs
  from database (fasta file):   'AT4G02541'     'AT1G12410.1'   'AT4G26510'
    'AT4G02540'     'AT3G04500.1'   'AT3G04500.2'   'AT3G04500.3'   'AT2G36250'
    'AT5G63350'     'AT2G36255'     'AT2G12280'     'AT1G03970.1'   'AT5G18390.1'
    'AT3G04490.4'   'AT5G38386'     'AT3G04490.2'   'AT3G04490.3'
    'AT3G04490.1'   'AT5G38380'     'AT4G02530'     'AT4G26520'     'AT1G69280.2'

link to the full stdout log after running the command without cds.fa
https://www.dropbox.com/s/gxa4e533p64wxm9/SnpEff_build-log.txt?dl=1


text:
Downloading SRA Files from AWS
I want to download the original BAM files that the authors had uploaded to SRA. Normally, I would just use sam-dump, but the files are having issues that seem related to this issue. Since according to the entry, AWS S3 also hosts the original BAM files, I thought I could download these directly.
Example SRA Entry

NCBI documentation implies that I can't download this directly, but I can freely copy to other AWS locations within the region. To this end, I created my own S3 bucket (mm-mneuron) and am now trying to copy from the SRA bucket to mine.  Here's what I try:
import boto3
import botocore

s3 = boto3.resource('s3')

bam_file = {
  'Bucket': 'sra-pub-src-6',
  'Key': 'SRR5253957/RPI25_0.bam'
}

my_bucket = s3.Bucket('mm-mneuron')

my_bucket.copy(bam_file, 'RPI25_0.bam')

This fails with:

botocore.exceptions.ClientError: An error occurred (403) when calling
  the HeadObject operation: Forbidden

That is, it sounds like I can't access the SRA bucket. I've tested downloading and uploading to my bucket, so I know I have write permissions. Not sure what else to try here.
How can I access the SRA data on S3?

It seems like is not a public bucket. You should download inside an (SRA AMI)[https://www.ncbi.nlm.nih.gov/sra/docs/sra-aws-download/] from sra-pub-src-6 to your instance, and then, in you local terminal download from SRA AMI in EC2 to your computer.

@zorbax yes that's what I did end up doing (see answer), but there is something specifically wrong with the deposits I'm trying to download. In general, are SRA buckets only accessible through SRA tools? i.e., one has to run an EC2 instance? It seems like direct transfers across S3 buckets would an obvious interface to provide.

You only have access to the public data sets with s3, and the key prj_phs710EA_test.ngc in EC2 only bring acces to 1000 Genomes data. https://www.biostars.org/p/393172/

A member of the SRA submission staff pointed out that using
prefetch --type all SRR5253957

will download the original files. In this case, it means running the above within an EC2 instance colocated with the S3 bucket (so, us-east-1) and having installed and configured SRA Toolkit to work from AWS (as per this documentation).
Unfortunately, the particular files I am concerned with are not currently accessible, but this should generically work in most situations.


text:
What is the meaning of these misaligned reads in a sequencing run?
I am analyzing some SARS-CoV-2 sequencing runs abd often find read alignments like the one in the image.
GACCCGGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
GACCCTGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
          GTTTCAC---ACTTAAAAACACAGTCTGTAC
          GGTTTAC---ACTTAAAAGCACAGTCTGTAC
CCACACTCTCCTAGCAC---ACTTAAAAACACAGTCTGTAC
CCACACTCTCCTAGCACCATACTTAAAAACACAGTCTGTAC

It basically says there are reads that share part of a sequence and part of another one. The first 4 reads are the same sequence and it's a sequence that you can effectively find in the official SARS-CoV-2 assembly but the other two reads are a mix of this one with another sequence (and both are the same mix) but with only a small (3bps) indel of difference.
The sequence of the last two reads does not appear anywhere in the official SARS-CoV-2 sequence.

What does it mean?
Are these sequencing errors?
If that is the case how is it that the same mistake is repeated with similar but not exactly the same data?
If these are real RNA sequences present in the probe, why are these not used anywhere in the final assembly and how did they arrive there?

Is there any other option that I am missing?

That's an old question that somehow has reappeared today. I'd like to point out that I realized months later that there was a problem in my alignment software and I think that's why there was a repeated indel.

There is a lot of biological interest in this sort of question, so its really important to understand this thoroughly

The unaligned parts of the bottom two sequences are from elsewhere in the sars-cov-2 genome if you blast them. Likely causes of this are library prep. artifacts (not terribly common, but you'll find a few cases here in there once you have millions of reads) where fragments end up ligating together or (probably more likely) rearrangements of the viral sequence. My guess would be that the latter is what occurred, given how rapidly viruses reproduce. Whether that arrangement can lead to viable viruses I don't know (that's too far out of my field).
Assemblies represent at best a snapshot of a single example of a species at a single moment in time. Actual individuals in the wild can vary quite substantially, especially when they don't need to be viable.

Thanks Devon. Isn't it curious that there was a rearrangement at this location and next an error in replication at two bp distance?

Not really, locations with issues tend to be near to another.

I have news about that. I have found that sometimes there is a merge between a read and a reverse-complemented read. That can only mean that the merge occurs when the sequence has been translated to DNA in the PCR phase. Therefore at least in some cases it doesn't seem to be a rearrangement of the virus.
An example read for SARS-CoV-2:
AGGTAATAGCACATCACTACGCAACTTTAGATACATTTCTTTATTTAACAAAAAGGTGCACAGCGCAGCTTCTTCAAAAGTACTAGTTTAAAAGACCAATAAATCCTACTGACCAGTCTTCTTACATCGTTGATAGTGTTACAGTGAAGAATGGTTCCATCCATCTTTACTTTGATAAAGCTGGTCAAAAGACTTATGAAAGACATTCTCTCTCTCATTTTGTTAACTTAGAC

That happens during library prep whenever the amplified fragment is less than twice the read length.


text:
TEs vertical or horizontal transfers?
After predicting and annotating all the transposable elements. For four related species.
Now I would like to find out these predicted candidates found in my costume library, either vertical or horizontal transfers?
are there any methodologies to suggest?

You make a phylogenetic tree of the "host species" and compare it to the phylogeny of the TE's. Homology is not a good idea for evolutionary trees.
There are two snags:

Concerted evolution, the TE's will risk recombining with each other within a host genome and that can mess up the phylogeneic signal (the theory as why is quite complicated)
If host species phylogeny are higher eukaryotes thats a lot of genetic material, so you'd need a short-cut. You can either use someones species tree, or focus on TE's within a subsection of the host genome and produce a tree of this subsection.

With point number 2 "a subsection of the genome" you can then assess synteny, i.e. has the TE jumped between position between species. If the TE hasn't moved then its an orthogonal proof of putative vertical transmission.
What you are looking for is incongruence between host phylogeny and TE phylogeny as evidence of horizontal transmission and congruence to call vertical transmission. However, incongruence could occur due to the complexity of concerted evolution, so synteny would assist again as an orthogonal proof. Congruence is assessed as a "mirror tree", there is quite a lot of theory on this looking at host-parasite coevolution
For making trees RAxML or PhyML are flavours of the month particularly for TEs.

A state of the art-of-the-art review is here, but it goes into significant detail into precise proofs, which in my opinion arn't needed

Thank you for the information, would you please provide a reference or study I can rely on. help me with the confusion

Sure @bioinfo, I'll do this at the beginning of next week. It will be pathogen coevolution rather than TEs, but its the same idea. Please remind me here if I forget

The review I would recommend is here https://academic.oup.com/cz/article/62/4/393/1745416 . However, the techniques are heavy duty state-of-the-art. I'll look at simpler approaches because you don't need this level of proof in my opinion


text:
Assembly reads having a copy of their beginning in their tail
I am analyzing the reads for the SARS-CoV-2 assembly with id SRR11140748. Apparently these reads were obtained with parallel sequencing by Illumina and Oxford Nanopore Technologies.
I have found these reads
GGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGGGTTTGCGGTGTAAGTGCAGCCCGTCTTACACCGTGCGGCACAGGCACTAGTACTGATGTCGTATACAGGGCTTTTGACATCAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGAGTATGTACAAATACCTACAACTTGTGCT

As you can see, the second one has an exact copy at the tail of its beginning:
AGTATGTACAAATACCTACAACTTGTGCT

I have found other issues apparently caused by rearrangements of the virus genome, but this case seems to be different. It seems unlikely that the virus mutated in such a way and then a read was cut exactly this way. What is the most likely explanation for this read?
It's not the only one displaying the same or similar behavior. for example this other one has a copy at the tail but from bps not exactly at the beginning (CACCGGAAGCCAATATGGATC):
CAGTTACACCGGAAGCCAATATGGATCAAGAATCCTTTGGTGGTGCATCGTGTTGTCTGTACTGCCGTTGCCACATAGATCATCCAAATCCTAAAGGATTTTGTGACTTAAAAGGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTTACACCGGAAGCCAATATGGATC

Did you check if that repeated sequence perhaps corresponds to a barcoding adapter?

I can't check that. It's not my data and I don't know what barcodes they used. But I guess if it was a barcode it would be repeated everywhere. I have made a BLAST search of some of these sequences (CACCGGAAGCCAATATGGATC for example) and they are part of the virus sequence.

If you google AGTATGTACAAATACCTACAACTTGTGCT you'll find it's a primer sequence: https://artic.network/resources/ncov/ncov-amplicon-v3.pdf

I still don't get it. Even if that is a primer I don't see how it explains it being repeated at the head and tail.

I guess it means the cDNA generated in the PCR sometimes merges with another cDNA.

Chimeric molecules are known to happen yes due to either PCR or the following ligation step.


text:
Is there a way to save a full pdb file in PYMOL after uploading biological unit?
pdb code: 2obt
The pdb file has a coordinates of a monomer.
I have used following commands to generate an assembly, or in other words to upload one more biological unit so that the protein would be a dimer:
set assembly, ""
fetch 2obt, asu, async=0

set assembly, "1"
fetch 2obt, assembly1, async=0

After generating the functional dimer in  pymol under name assembly1, I have an important issue that I cannot save a dimer structure in pdb file with dimer coordinates. Does anybody has an idea how to go around it or is it even possible to do that?

can we says that "after uploading the biological unit is not the right wording" ?? --> When loading mmCIF structure files, the assembly setting controls whether PyMOL loads the asymmetric unit (assembly="") or a biological unit (e.g. assembly="1"). Assembly identifiers in mmCIF files can be arbitrary strings, but usually are numeric and most files define at least assembly "1".  I guess the biological unit is created applying same symmetry operator on the asymmetric unit (or am I wrong ??). Ok PymolWiki says loading, actually loads .....

So on the GUI form of PyMOL on the right hand side you have a list of what is ambiguously called objects in PyMOL (models elsewhere), which can be enabled or disabled by clicking on their name.
The name also has a 1/1 in green or similar. This indicates the state. In your case you have 2 states.
To split a multistate object use
split_state prot_A

To combine two objects one uses the command
create new_name, prot_A or prot_B

Now, this has two caveats which if violated result in a nasty piece of modern art:

if the chains, segment ids and residues are the same
if you had two overlapping macromolecules

To deal with the first, change the chain id on one of them.
alter prot_B, chain='B'
   sort
Segment identifier is a weird thing that most comp-biochem tools dislike, so I like to remove that
alter prot_B, segi=''
   sort
Now when you go to save your new_name object, in the popup of the GUI you get a dropdown menu defaulted to enabled (visible objects that have not be clicked to be disabled). There you can choose new_name.
Aligning two
If someone that is in not the OP author has a template protein and needs to make a dimer of another based on that, the process is the same.
create twin, prot_B
align twin, prot_A and chain B
alter twin, chain='B'
sort

SDF files
If someone else has a multicompound sdf where the state is the header name, I strongly advise not to use split_state. Yes, the python pymol module (not the app) handles states in a really weird way, but states can be selected with state 1 selection algebra and state is an argument for most commands that do not accept state selector.


text:
conda doesn't install latest version of snakemake
I would like to use the --default-resources parameter for profiles which is available in later versions of snakemake.
To install snakemake, I created a new conda environment:
conda create -n snakemake -c conda-forge -c bioconda snakemake
The currently installed version is listed as 5.14:
$ conda activate snakemake
(snakemake)$ conda list | grep snakemake
# packages in environment at /home/uvi/be/hde/.conda/envs/snakemake:
snakemake                 5.14.0                        0    bioconda
snakemake-minimal         5.14.0                     py_1    bioconda

However, command line tool says it is on version 5.1.4:
(snakemake)$ snakemake --version
5.1.4

This version does not know the new parameter --default-resources yet. So it gives an error when using it

snakemake: error: unrecognized arguments: --default-resources partition = "thin-shared"

I already checked that I'm running the executable from inside the environment, rather than a different, external one:
(snakemake)$ which snakemake
~/.conda/envs/snakemake/bin/snakemake
(snakemake)$ $(which snakemake) --version
5.1.4

What is the output of python -c "import snakemake; print(snakemake.__file__)"?

/home/uvi/be/hde/.local/lib/python3.6/site-packages/snakemake/__init__.py Hah! That's not the version inside the conda environment.

Python (whether installed with conda or not) will prefer packages installed in your home directory unless you tell it not to. This is a little known secret about using a virtualenv, actually, since internally those tell python to ignore everything in your home directory. Within conda, there's long been a debate about whether python (and also R) should continue looking for packages in the normal places or only within a conda environment. Since it would likely break many people's workflows, it's been decided to allow the normal language features of looking for packages not just within conda envs but also in the normal system and user locations. For python in particular you can disable this as follows:
PYTHONNOUSERSITE=True snakemake ...

Then only the packages installed within your conda environment will be used.

This is very useful, thanks a lot @DevonRyan! I solved the problem by uninstalling the external package (which I found thanks to your comment above). But for the future I will consider using PYTHONNOUSERSITE=True.

https://xkcd.com/1987/


text:
Extract end position from CIGAR
I was wondering how could I obtain the end position of an alignment using its start position and its CIGAR String.
If an alignment present some soft clipping, for example:

    20S55M

Then, would the start position indicated in the SAM file be the first soft clipped base, or the first correctly mapped base?
Thanks

The start position is the first correctly mapped base. One way to do it, is to use GenomicAlignments in R to read in the alignments.
For example the sam file looks like this:
cat test.sam
@HD VN:1.3  SO:coordinate
@SQ SN:chrI LN:228861
D00689:415:CD4LUANXX:3:1316:7447:98083#CGGCTATGCAGGACGT 0   chrI    769 0   42M8S   *   0   0   CCATAGCCCGGTGGGGCCTGTGCTGCCACAGAGGGTGTGTCTAGATCGGA  B/B//FFFFFFBFFFF//7FF/<F/FFBFFFF//<F///<F/BFB//<F/  NM:i:1  MD:Z:21A20  AS:i:37 XS:i:37
D00689:415:CD4LUANXX:3:2308:9597:30445#CGGCTATGCAGGACGT 0   chrI    3965    0   3S47M   *   0   0   CCACACACCCACACACCACACCCACACCACACCCACACCACACCCACACC  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFBF<FFFFF<FF  NM:i:0  MD:Z:47 AS:i:47 XS:i:47

We convert to bam:
samtools view -bS test.sam > test.bam

And read using R:
library(GenomicAlignments)
gr = readGAlignments("test.bam")
gr

GAlignments object with 2 alignments and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chrI      +       42M8S        50       769       810        42
  [2]     chrI      +       3S47M        50      3965      4011        47
          njunc
      <integer>
  [1]         0
  [2]         0
  -------
  seqinfo: 1 sequence from an unspecified genome

The start and the end can be retrieved using:
start(gr)
[1]  769 3965
end(gr)
[1]  810 4011

The match is starting with the first mapped base ignoring any softclipped parts of the read.
As an example this are two reads from a bam file of mine:
start   cigar   sequence
11710   17S83M  NNNNNNNNNNNNNNNNNGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCC
11713   100M    CTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCT


text:
Chromosomal order of supplementary alignment in BAM file "SA" tag
I have several long-read sequences, which when aligned to a human reference genome, aligns to multiple chromosomal fragment as a result of chromothripsis shattering the linear DNA and randomly piecing them together. This creates for a single primary alignment multiple supplementary alignments.
With a single example of a read (only showing parts of all the bamfile information) being
cdaf1fcb-27bf-4890-8ac7-a097dfa12c38    16      chr11   15151487        60
With the primary alignment starting at chr11 at 15151487
Then the supplementary alignment can be seen as follows
SA:Z:chr13,49827969,+,1839S956M25D1584S,60,53;chr11,15152459,-,833S749M22D2797S,60,72;chrX,138324656,-,2540S404M2D1435S,60,20;chr11,15152066,+,278M11D4101S,60,20;
Where the alignments are to chr13, 11, X and 11. With both regions on chr11 being close to the initial region from the primary alignment, making the regions on chromosome 11 the same (I know I haven't calculated the full chromosomal interval using the Cigar string for this example).
I assume the order of that supplementary information above is in the actual order. But the problem is that my DNA should be circular which doesn't fit with the order of the alignments since the circle would then be chr 11 (primary, start of circle) -> chr 13 -> chr11 (same as primary) -> chrX -> chr11 (back to the start of circle)
My question is:
Q) So is the order of the supplementary information just pieces together randomly or is it true that I can assume the order is the true order of the sequencing read?
So thanks a lot for your help.

There is nothing in the SAMtags specification that dictates the order of the entries in an SA tag.  It may be entirely dependent on the aligner that you used.
Looking at my own data that was aligned with bwa, I cannot see a definite pattern.  Longer matches tend to be first, but not always.
Your aligner documentation may have more information on how it uses SA tags.


text:
General question: How to assign my SNP (Wheat, Illumina 15K) to specific chromosomes?
I have 2000 lines of wheat genotyped with a 15K illumina chip.
I also have the name of the SNPs.
I would like to know how could I find out which SNPs are in each chromosomes.

Answer from @devon-ryan, converted from comments:
Look in Biomart to see if there's SNP rsids for wheat.


text:
Significance of k-mer length in COVID-19 sequence analysis?
I'm getting started in biology and bioinformatics with sequencing the SARS-Cov-2 or Coronavirus genome.  I'm interested in this code which identifies k-mers in the genome:
from skbio import DNA
DNA(str(seq.seq)).kmer_frequencies(5)
kmers = DNA(str(seq1.seq)).kmer_frequencies(5)
kmers
{'ATTAA': 60, 'TTAAA': 95, 'TAAAG': 64, 'AAAGG': 45, 'AAGGT': 46, 'AGGTT': 51, 'GGTTT': 55,
'GTTTA': 69, 'TTTAT': 75,...}

The author of the code wants to use k-mer breakdowns to vectorize the data for machine learning purposes.  What I'm curious about is why the k-mer in this case is five nucleotides?  What factors should go into my selection of the length of the k-mer for analysis?

5bp produces a high resolution analysis, providing its sufficiently robust.
Its a sliding window of 5 bp wide and 1 bp step length which runs along the length of each genome. They are performing SNP analysis/classification. The original sequence they were examining was,
ATTAAAGGTTTAT

In any sliding window analysis the window size is arbitary towards the resolution and robustness of the analysis. The author believes that 5bp provides sufficiently robust information for a very high resolution classification/training analysis in their example.
The proof is the accuracy of the classification. I assume they have good results from 5bp. In ML the rationale is the general laws of ML (over/under training) and assessed on the test split in the data (probably 20%).
A key reason for the 1bp step is that this is the way ML works, you need a lot of data and by doing 5 bp sliding window at 1 bp step that will generate ~30 000 pieces of data for each sequence. So that's 3 million data points for 100 sequences.

Would the SNP analysis you mention include a classification algorithm such as principal component analysis?

Oh sorry thats different, PCA (ML) is a different analysis to the one I had in mind (supervised learning) and is usually performed in conjunction with tSNE @greenpenguin. PCA is judged on its ability to cluster the classification targets into different groups. PCA/tSNE is called unsupervised learning


text:
Programmatically retrieve graph features from Cytoscape
Is there a way to programmatically retrieve graph features (a networkx graph) from Cytoscape in Python 3?
EDIT: To add more details, I had a dictionary of proteins, each protein was associated with a list of proteins that it interacts with. I wanted to create a protein-protein interaction graph so that I could extract graph features from it, including things like k-core value of each node.

Welcome to Bioinformatics Beta! It would be helpful if you added more details. What have you attemped? What did you observe when you made the attempt that prevented you from getting the answer you wanted?

Extract this information into what? A Python dictionary? JSON? Pandas DataFrame? It might help the community understand what your issue was if you elaborate.

I was able to extract the graph features that I wanted just using networkx, without going through the intermediate step of Cytoscape.
EDIT: Initially, I planned to import the data into Cytoscape and extract graph features from Cytoscape. However, I noticed that a much easier solution would be to just use the methods of the networkx library. I initially had a few graph features in mind that I wanted to extract like k-core value of each node, and I was able to do this using methods from networkX. To calculate k-core for instance, I used the k_core method. For more information on networkX methods, check the docs

I'm glad to hear that you solved your own question. If you edit your question to clarify what the issue was, and edit this answer on how you solved it, then this question can provide value to others who come across the same issue.

Implementation details would be useful. For example, did you use nx.get_node_attributes, nx.get_edge_attributes, or something else?


text:
Protein misfolding
I am looking for literature on protein and RNA mis-folding. I am on the computational/biophysics side, and what interests me are the practical applications: diseases caused by misfolding, drug discovery, etc.

Protein mis-folding is a very interesting area right now @vadim

https://en.wikipedia.org/wiki/Elexacaftor/ivacaftor/tezacaftor

Excellent question! I would recommend that you start with an NIH resource at the National Library of Medicine (NLM). PubMed (https://pubmed.ncbi.nlm.nih.gov/) provides search capabilities across the biomedical literature. You can also create an account and use it to manage your searches of the literature. An example search below for "protein misfolding disease" in PubMed retrieved almost 50,000 hits. You can filter to make the results more specific to your interests. You can view the citations, including the abstracts, and often links to the full text. Since you mention that you are on the "computational/biophysics side," you might be interested in the PubMed API (https://www.ncbi.nlm.nih.gov/pmc/tools/developers/) and python modules that use the API, such as Biopython.
E-utilities:

https://www.ncbi.nlm.nih.gov/books/NBK25501/

Biopython:

https://biopython.org/DIST/docs/api/Bio.Entrez-module.html

PubMed Search:

https://pubmed.ncbi.nlm.nih.gov/?term=protein%20misfolding%20disease

Specific Articles:

https://www.nature.com/scitable/topicpage/protein-misfolding-and-degenerative-diseases-14434929/
https://www.nature.com/articles/nrn1007


text:
How to generate the json files for Cromwell workflow execution?
I am starting to use Cromwell WDL scripts for running GATK workflows on a slurm cluster. When using the paired-fastq-to-unmapped-bam script I need to define the fastq sepcifics in the appropriate json file. This needs to be written separately for each pair of fastqs. My question is: Is there a way of generating the json files on a batch of fastq-files by scripting? Of course, I could start writing scripts myself, but since I assume that GATK is used in a batch environment, I have the hope that something already exists.

Hi @luckyman, if you post the non-batch code you have written I'm pretty sure someone will loop it for you.

Unfortunately generating inputs for Cromwell is not something that is handled by Cromwell itself or by GATK, and generally needs to be scripted to fit your custom needs.
Despite this being a bit arduous, there is a very good reason for this. WDL itself is a highly polymorphic workflow language, where inputs can be vastly different between each workflow. Cromwell ONLY understands how to parse the workflow and the inputs, but does not know the context of what it is running (ie this is a paired mapping job), nor does GATK know that it is running within a WDL.
Additionally, I've encountered many different cases where there are not just 2 fastq files but, due to older sequencing technologies, 10 or more for the same sample. Unless GATK/Cromwell knows how your specific samples are arranged there is no easy way for them to generate your inputs.
There are tools that would possibly help you simplify scripting better. WOMTool developed by the Broad has an inputs command which will take an  arbitrary workflow and generate a JSON template of the inputs, which you can then use in your own scripts.
As an example, imagine you have this WDL workflow in a file called my.wdl:
version 1.0
workflow foo {
  input {
    Array[Pair[File,File]] fastqFiles
  }
}

When you run WOMtool's inputs command, it produces a template JSON listing all inputs and specifying the variable types expected for each:
$ java -jar womtool.jar inputs my.wdl
{
  "foo.fastqFiles": "Array[Pair[File, File]]"
}

You can use this template to produce actual JSONs for each of your samples.
Note that you can also get the template written to a file as follows:
$ java -jar womtool.jar inputs my.wdl > my.inputs.json


text:
How to calculate average minor allele frequency (MAF) difference between two populations (i.e. European, Hispanic)
The data is present in the PLINK format .bed (individual genotypes), .bim (genetic markers) and .fam (sample IDs and disease phenotype file). I have also separate text files containing information on subjects related to two populations. (i.e. European and Hispanic)

You need to report allele frequencies using plink for some European samples
Do the same for Hispanic samples
Use a software to subtract the two allele frequencies obtained above and obtain the result.

Hi, I think I kind of understood what you mean. So in PLINK, we use --freq command to find information on frequency. and what software do you recommend to subtract the frequency files of two populations? thank You


text:
Generate ligands candidates based on protein shape
Recent approaches to novel drug design using machine learning (ML) and deep learning, often involve generating hundreds of potential ligands which are later tested by docking with a target protein and recording the resulting binding affinity.
Is it possible to go the other way around by first constraining the ligand generation process to the shape of the protein thus limiting the search space, as opposed to generating ligands of random shapes?

For an example of an approach that builds up fragments in pockets to find ligands, see FTMap: https://www.nature.com/articles/nprot.2015.043

Thank you @jgreener, this is the kind of approach I was looking for!

TL;DR: docking is much slower than any ML approach, but the ML approach can be constrained by pharmacophores dictated by the active site.
Side note: Scale
The scale for ligand space exploration is generally several orders of magnitude higher than "hundreds": Zinc DB lists 750 million enumerated compounds, GDB13 has enumerate all possible compounds up to 13 atoms (970 million), Enamine Real is 1.2 billion etc. In an actual experiment, one would use a more restrictive set. As an example of a real project, in Covid19 moonshot project there are 2,800 user suggested follow-ups.
Side note: More than just binding
There are a few requirements for a drug to be effective, and binding is only one of them:

Must bind...
Must follow Lipinksi's rule of five, i.e. it must be membrane permeable and small.
Must be synthesisable. Whereas the experience of someone learnèd is by far the best way to tell how an expansion can be done, there are several machine leaning approaches that can filter a huge dataset. IBM actually have a fun online tool —.
Must not have side effects

So pre-filtering the list of compounds, regardless of docked poses, will happened, be it manually, by cut-offs or by advanced machine learning.
CPU time: docking vs. machine learning
These various steps are generally followed in that order.
Calculating properties (MW, logP etc.) for a given SMILES string to filter based on Lipinski's rule takes a fraction of a second, so that is generally always a first step.
Whether the compound can synthesised or bought is a general left last and done manually. Pharmacology is complicated so is ignored.
Docking a ligand takes time and different approaches have different time scales. A decent algorithm that uses implicit solvent model takes about a minute per core. A short MD simulation takes an 1 hour per core or more. There are docking programs for screening that can be faster but these dock only one conformer of the small molecule against a rigid protein and the results are not of much value. Also setting up docking properly is not easy (cf. steps required for Autodock 4, each with their own caveats). So to recap what needs to be considered that results in the slowness:

Hundreds of conformers of the ligand
Some flexibility of the protein, side-chains repacking only or even backbone changes.
thousands of poses explored.
Optionally. Implicit water is not as good as a TIP3 water model

Machine learning have a much smaller time requirement and is way more attractive for publications. Generally these work by using ligand only information, such as properties, splitting it into pharmacophores, colour etc. This requires a decent dataset, such as empirical results, say from a fragment screen, what compounds bound and which didn't. Often there is a lot of hype with machine learning and the result is something close to something already seen, as well put in this blog post by Pat Walters. This is a side effect of the starting data. The recent Stokes et al. 2020 Cell paper uses AI perfectly because of their staggering automated empirical data collection pipeline.
ML after docking
However, ML can actually be used after docking to improve the quality of the score.
The Boyles, Deane, Morris 2019 Bioinformatics paper investigated ligand-only machine learning to see why do ML methods that rely solely on ligand somehow often work. It is a good paper from a sceptical computational biochemistry viewpoint.
The scorefunction used in docking by the sampler can be classified as forcefield-based, empirical or hybrid. In some cases empirical cases, many terms are ligand based properties (/features) that would normally go into a ML search (say one trained on failed hits and successful hits).
For example Autodock Vina (not 4) has regression derived factors in the final weights of the score
Structural features in machine learning: HotSpots
Information of the protein active site can be used in ranking compounds in a non-docking machine learning approach. Namely, finding out what kind of chemical groups bind preferentially to different parts of the protein's active site, i.e. making a "fragment hotspot map". Specifically, a hotspot map allows the determination of a pharmacophore, which is an abstract generalisation of a group of molecules with both both electronic and steric features —a coarse-grain model, if you will. This approach is also used to detect where the active site actually is in the extreme cases where it is not know (primarily a problem for pipelines, not humans). There are several implementations as described here, the first being FTMap (site, paper), which docks 16 different types of compounds. Waters can play a role and severely complicate matters, hence approaches can be taken to counter that. However, unlike docking protein rigidity is not a big nuisance: this is because a hotspot map is most often used for filtering/weighing by the "color" not the shape of the molecules. But shape based searching is possible.
Shape based scoring
In terms of shape-based searching, this is actually a question of having pre-computed conformers, which can be problematic (but can and is often done on the fly), and shape based scoring —for which there are different metrics, such as the classical RMSD to more complexly weighted scores such as SuCOS.

A very comprehensive and detailed answer @matteroferla

Your answer seems to reinforce the idea that docking ligands of random shapes with target proteins is a slow and often complicated procedure, if so, wouldn't it make more sense to first identify binding hot spots and take the surface geometry of the target protein into account before generating any ligands. It seems to me that this approach would greatly limit the search space by filtering out ligands that don't match/fit with the surface geometry of the protein? Thank you!

Busted! I had actually considered talking about using HotSpots to filter out pharmacophores of interest, but I did not want to explain pharmacophores and it is not a too mainstream method. There are few implementations, but Radoux and Blundell, 2015 is a good read. I'll update the answer when I get a chance.

Thank you @MatteoFerla, I really appreciate your input, I have set your answer as accepted.

Thanks, I have amended a bit my answer, but I feel like I kind of missed your question in the first instance —I thought it was more regression/machine learning focused...


text:
What's the meaning of having different alignments with the same qname in a bam with no duplicates?
I have removed duplicates from my bam file (not secondary, nor suplementary alignments), and then I have looked for some repeated qnames (to assess the failure of a program) and I've found a lot of them. All of them has TLEN equal to 0.
As an example: this is part of a bed file that a generate from the bam.
CHR  Start   End       QNAME                                 PNEXT  TLEN
12 86945793 86945944 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0
4 168461445 168461485 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0

And those are these reads in the original bam:
A05325:240:SS348DSZX:2:5105:12246:3790  2177    4       168461445       0       111H40M 11      90044835        0       ATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA        :,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF     NM:i:2  MD:Z:28T8G2     MC:Z:151M       AS:i:32 XS:i:31 RG:Z:HM327DSXX.3        SA:Z:12,86945793,+,20S37M94S,1,1;
A05325:240:SS348DSZX:2:5105:12246:3790  129     12      86945793        1       20S37M94S       11      90044835        0       ATGGTGAAAAAACAAATTATTTAAAAAAATCAAAAAAAAAAAAAAAACCAAAAAACATAAAAAAAAGGGCAGTAATTTTATTCATGGGGTCAAACAAAAAAAAAGCATAATATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA      FF,,,:,F,F,,,,F,,,,,,,F,FFFF:,:FFFF:FFF,FF::F,:,,FFFF,,,F,FF::,,F,::,,,,F,,F:FF,,,,,,,F::F,F,,,,F,FFFFF:,F,,F,,:,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF      NM:i:1  MD:Z:9G27       MC:Z:151M       AS:i:32 XS:i:29 RG:Z:HM327DSXX.3        SA:Z:4,168461445,+,111S40M,0,2;

This are the commands used to remove duplicates:
samtools view -F 1024 -bo filtered.bam original.bam

The commands used to generate this bed were:
samtools view sample.bam | awk '{ print $3, $4, $4+length($10), $1, $8, $9}' > sample.bed

Why do I present those alignments?

Can you tell us more about how you generated the BED? I'm not sure that I know what TSLEN is off the top of my head. Googling around I found https://www.biostars.org/p/78878/ , maybe you are talking about TLEN (e.g. "TLEN: signed observed Template LENgth")? Have you already filtered alternate mappings from the BAM?

Those reads are from different tiles (3105 and 5105).

You're right! I've just edited to clarify these issues.

Please post the original SAM entries (before any filtering). Also describe exactly how deplicate removal was performed (with the commands).

Okay, i've added this information

I think that you are dealing with two different phenomena:
1) TLEN == 0 reads
From the BAM spec:

TLEN: signed observed Template LENgth. If all segments are mapped to the same reference sequence, the absolute value of TLEN equals the distance between the mapped end of the template and the mapped start of the template, inclusively (i.e., end − start + 1).14 Note that mapped base is defined to be one that aligns to the reference as described by CIGAR, hence excludes soft-clipped bases. The TLEN field is positive for the leftmost segment of the template, negative for the rightmost, and the sign for any middle segment is undefined. If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs. It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).

Note that in each case the mate read appears to be mapped to a different chromosome/contig. So, the mapper can't estimate insert size in these cases, so it just says zero.
2. Multiple reads with same qname
There are various reasons why the same qname can appear multiple times. Again from the BAM spec, one such explanation is that it is a chimeric alignment:

Chimeric alignment: An alignment of a read that cannot be represented as a linear alignment. A chimeric alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the others are called “supplementary” and are distinguished by the supplementary alignment flag. All the SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80 flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

You may have to do some digging and looking at the different flags to decide what exactly is going on, but my bet is that they are supplementary alignments.


text:
Intersecting a data frame based on the multiple columns
I have 2 data frames
> head(df1)
   Chromosome Start_Position End_Position       Gene
1:      chr12       46254706     46254706       ARID2
2:      chr17        7579849      7579849        TP53
3:       chr9       21974777     21974780      CDKN2A
4:       chr1       27105930     27105930      ARID1A
5:       chr1      201982973    201982973        ELF3
6:       chr3       89480459     89480459       EPHA3
> head(df2)
  Chromosome Start_Position End_Position
1          1          10583       863511
2          1       12841835     12854479
3          1       17231299     17232877
4          1      120380739    142540462
5          1      142540484    142684781
6          1      142685823    143462598
>

How I can intersect genes from df1 which are in the range of df2? Because both data frames have chromosome number, start and end but I don't know how to identify genes from df1 based on df2
I have tried
df1[df1$Start_Position%in%df2$Start_Position,]

But I obtain one gene
I also tried
> new_dataset <- df1 %>% right_join(df2, by=c("Chromosome","Start_Position"))
Error: Can't join on 'Start_Position' x 'Start_Position' because of incompatible types (character / integer)
>

reformat df2 so that you have chr1 is instead of 1, then you basically have two bed files.  Then use bedtools intersect.  https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

In R:
> write.table(df1, 'df1.txt', quote=F, sep="\t", row.names=F, col.names=F)
> write.table(df2, 'df2.txt', quote=F, sep="\t", row.names=F, col.names=F)

Convert both files to sorted BED:
$ sort-bed df1.txt > df1.bed
$ awk -v FS="\t" -v OFS="\t" '{ print "chr"$1 }' df2.txt | sort-bed - > df2.bed

Then map df1 gene names to df2 coordinates:
$ bedmap --echo --echo-map-id-uniq --delim '\t' df2.bed df1.bed > answer.bed

To bring back into R, e.g.:
> dfMap <- read.table('answer.bed', sep='\t', col.names=c('Chromosome', 'Start_Position', 'End_Position', 'Genes'), colClasses=c("factor", "numeric", "numeric", "character"))


text:
Grouping a set of genes in a ggplot object
I have a boolean matrix of mutated or non mutated genes across some samples
> head(muts)
                  TP53 CDKN2A TSHZ3 ARID1A EPHA3 NOTCH1 ACVR2A ARID1B ARID2
LP6008336-DNA_G01    1      1     1      0     0      0      0      0     0
LP6008460-DNA_D01    1      1     0      1     0      0      0      0     0
LP6005500-DNA_D03    1      1     0      0     0      0      0      0     0
LP6008334-DNA_B02    1      0     1      0     0      0      0      0     0
LP6008334-DNA_A04    1      0     1      0     0      0      0      0     0
LP6005334-DNA_H01    1      0     0      0     1      0      0      1     0

By this code
pw_rt = data.frame(
  gene = c('PIK3CA', 'EGFR', 'ERBB2', 'PTEN'),
  pw = c('RTK'),
  stringsAsFactors = F
)

pw_cc = data.frame(
  gene = c('CDKN2A', 'CCNE1'),
  pw = c('Cell Cycle'),
  stringsAsFactors = F
)

pw = rbind(pw_rt, pw_cc)

pw_col = RColorBrewer::brewer.pal(n = 3, name = 'Set1')[1:2]
names(pw_col) = unique(pw$pw)

require(dplyr)

lmuts = reshape2::melt(muts) %>% as_tibble()
colnames(lmuts) = c('sample', 'gene', 'value')

pw_n = pw$pw
names(pw_n) = pw$gene

Np = nrow(muts)

lmuts %>%
  mutate(gene = paste(gene), PW = pw_n[gene]) %>%
  group_by(gene, PW) %>%
  summarise(N = sum(value == 1)) %>%
  ungroup() %>%
  mutate(
    PW = ifelse(is.na(PW), "None", PW),
    N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
    ) %>%
  ggplot(aes(x = PW, y = gene, fill = PW)) +
  geom_tile() +
  geom_text(aes(label = N)) +
  theme_light() +
  theme(legend.position = 'bottom') +
  guides(fill = guide_legend('Pathway')) +
  scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
  scale_x_discrete(limits = c(names(pw_col), 'None')) +
  labs(
    x = "",
    y = 'Gene',
    title = "Occurrence of mutations"
  )

I am able to get proportion of altered pathway by these genes like

However I want to highlight oncogenes and tumour suppressors within this plot
For instance TP53 is a tumour suppressor and CCNE1 is an oncogene
I don't know how to manipulate the above code to highlight or partition genes to tumour suppressor or oncogene
Can you help?
EDITED
I have add extra column to lmuts like
> head(lmuts)
             sample gene value gene_class
1 LP6005690-DNA_H02 TP53     1        TSG
2 LP2000333-DNA_A01 TP53     1        TSG
3 LP6005409-DNA_D03 TP53     1        TSG
4 LP6008141-DNA_H02 TP53     1        TSG
5 LP6008336-DNA_E02 TP53     1        TSG
6 LP6008269-DNA_B06 TP53     1        TSG
>

but
> lmuts %>%
+     mutate(gene = paste(gene), PW = pw_n[gene]) %>%
+     group_by(gene, PW) %>%
+     summarise(N = sum(value == 1)) %>%
+     ungroup() %>%
+     mutate(
+         PW = ifelse(is.na(PW), "None", PW),
+         N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
+     ) %>%
+     ggplot(aes(x = PW, y = gene, fill = PW)) +
+     geom_tile() +
+     geom_text(aes(label = N)) +
+     theme_light() +
+     theme(legend.position = 'bottom') +
+     guides(fill = guide_legend('Pathway')) +
+     scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
+     scale_x_discrete(limits = c(names(pw_col), 'None')) +
+     labs(
+         x = "",
+         y = 'Gene',
+         title = "Occurrence of mutations"
+     )+facet_wrap(~gene_class, ncol=1)
Error: At least one layer must contain all faceting variables: `gene_class`.
* Plot is missing `gene_class`
* Layer 1 is missing `gene_class`
* Layer 2 is missing `gene_class`
Run `rlang::last_error()` to see where the error occurred.

If you had the classification into oncogene/TSG for each gene in a column (let's call it gene_class) in lmuts, you could split both groups into facets within the plot by adding:
+ facet_wrap(~gene_class, ncol=1)

I can't reproduce your graph, because you don't make muts available, but it appears your problem is having x = PW in your ggplot call. Don't you want to that to be whatever your variable name for sample is (ie, the rownames in your but dataframe you have at the beginning).

Sorry I have edited my post, know you can reproduce my plot because muts is available now


text:
A .gtf file for HG38 for HISAT2 indexing
I am trying to build a HISAT2 index using Pertea et al. procedure (https://www.ncbi.nlm.nih.gov/pubmed/27560171) using their extract_splice_sites.py and extract_exons.py pieces of code. For that, I need a .gtf file of Human Genome HG38. I have obtained this file using recommendations described here http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download
The file I have downloaded is GCF_000001405.39_GRCh38.p13_genomic.gff (1,1 Gb) and when I use on it this command python3 extract_splice_sites.py GCF_000001405.39_GRCh38.p13_genomic.gff > test.ss my test.ss file is 0 Kb! I wonder if there's another .gff file that I must use and where I can get it? This is my first time working with the human genome so I am not sure if I am doing the procedure the right way.

you might wanna try python2.7 https://www.biostars.org/p/293881/

You can download a GTF file from NCBI now. For hg38 (also known as GRCh38), just go to NCBI Assembly and search for 'GRCh38'. Click on the first result to go to this page and use the blue 'Download Assemblies' button to download data. Here, after you choose 'RefSeq' as your source of data, you will see a GTF file to download from the drop-down list, as shown in the screenshot below:

I had the same problem and I was able to solve it with that
https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff_to_gtf.md