text:
How do I include repeat purity, default slippage, default stutter, and minimum flanking (left and right) in Tandem Repeat Finder's output?
I am attempting to create a markerInfoFile for use in a program called popSTR (GitHub Documentation: https://github.com/DecodeGenetics/popSTR). The marker info file should contain information about microsatellites in a reference genome. Microsatellites are stretches of DNA where a short motif (i.e. "AGG") is repeated multiple times (i.e. three repeats of "AGG" is "AGGAGGAGG"). Here is an example line from the popSTR author's marker info file:
chr10 10589 10601 GCCC 3.2 NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCACCCTTCTAACTGGACTCTGACCCTGATTGTTGAGGGCTGCAAAGAGGAAGAATTTTATTTACCGTCGCTGTGGCCCCGAGTTGTCCCAAAGCGAGGTAATGCCCGCAAGGTCTGTGCTGATCAGGACGCAGCTCTGCCTTCGGGGTGCCCCTGGACT GGTCTGTGCTGAGGAGAACGCTGCTCCGCCTCCGCGGTACTCCGGACATATGTGCAGAGAAGAACGCAGCTGCGCCCTCGCCATGCTCTGCGAGTCTCTGCTGATGAGAACACAGCTTCACTTTCGCAAAGGCGCAGCGCCGGCGCAGGCGCGGAGGGGCGCGCAGCGCCGGCGCAGGCGCGGAGGGGCGCGCCCGAACCCGAACCCTAATGCCGTCATAAGAGCCCTAGGGAGACCTTAGGGAACAAGCATTAAACTGACACTCGATTCTGTAGCCGGCTCTGCCAAGAGACATGGCGTTGCGGTGATATGAGGGCAGGGGTCATGGAAGAAAGCCTTCTGGTTTTAGACCCACAGGAAGATCTGTGACGCGCTCTTGGGTAGAGCACACGTTGCTGGGCGTGCGCTTGAAAAGAGCCTAAGAAGAGGGGGCGTCTGGAAGGAACCGCAACGCCAAGGGAGGGTGTCCAGCCTTCCCGCTTCAACACCTGGACACATTCTGGAAAGTTTCCTAAGAAAGCCAGAAAAATAATTTAAAAAAAAATCCAGAGGCCAGACGGGCTAATGGGGCTTTACTGCGACTATCTGGCTTAATCCTCCAAACAACCTTGCCATACCAGCCCATCAGTCCTCTGAGACAGGTGAAGAACCTGAGGTCGCAGGAGGACACCCAGAAGGTCCAGAGAGAGCCTCCTAGGCCCCCCACCTCCCCCCGTGGCAGCTCCAACCCCAGCTTTTTCACTAGTAAGGCAGTCGGGCCCCTGGGCCACGCCCACTCCCCCAAGCGGGGAAGGAGCTTCGCGCTGCCGCTTGGCTGGGGACTGGGCACCGCCCTCCCGCGGCTCCTGAGCCGGCTGCCACCAGGGGGCGCGCCAGCGGTGTCCGGGAGCCTAGCGGCGCGTGTGCAGCGGCCAGTGCACCTGCTCTGGCCCTCGCCGCGGTCTCTGCCAGGACCCCGACGCCCAGCCTGACCCTGCCATTCAGCGGGGCTGCGGCTCCA GCCCGCCCGCCCG 4 4 1.00 0.0149855 0.949846 0 0.75 0.25 0

Though it may be difficult to tell from the above example, the desired output is a 17 column space-delimited file formatted like so:
  chrom startCoordinate endCoordinate repeatMotif numOfRepeatsInRef 1000refBasesBeforeStart 1000refBasesAfterEnd repeatSeqFromRef minFlankLeft minFlankRight repeatPurity defaultSlippage defaultStutter fractionAinMotif fractionCinMotif fractionGinMotif fractionTinMotif

Since I am working with non-human genomes, I cannot simply use the author's sample marker info files. I could not find explicit instructions for the creation of this markerInfoFile in popSTR's full journal article (https://academic.oup.com/bioinformatics/article/33/24/4041/2525679), so I emailed the authors. they instructed me to use Tandem Repeat Finder (in the Linux command prompt) with the following options:
./trf409.linux64 reference_genome.fasta 2 7 7 80 10 22 1000 -d -h -ngs > reference_genome_TRFmap

Here is the documentation for TRF:
https://blaxter-lab-documentation.readthedocs.io/en/latest/trf.html
I applied TRF to my reference genome, but the output was in the following format:
Chrom
    startCoordinate endCoordinate PeriodSize    numOfRepeatsInRef   SizeofConsensusPattern  PercentofMatchesBetweenAdjacentCopies   PercentofIndelsBetweenAdjacentCopies OverallAlignmentScore  fractionAinMotif fractionCinMotif fractionGinMotif fractionTinMotif Entropy repeatMotif repeatSeqFromRef    1000refBasesBefore  1000refBasesAfter

Here is an example of the first two lines of my current output file:
@ScyDAA6_1;HRSCAF=23
108 120 5 2.6 5 100 0 26 0 0 15 84 0.62 TTTGT TTTGTTTTGTTTT CAGTAAAGTCTTTCTTTCCTCTAACATAGAAAGTACTACTAGATTAGTGC TCTGTGTATGCTCTCTATTCTCAACCTCCAGATGCCCGTTCACACTGAGC
414 424 5 2.2 5 100 0 22 63 0 18 18 1.31 ATGAA ATGAAATGAAA TGAACTGGAATGTACATGATTGAATTTAAATTACTTCTTTAAAAAATTCC AGTATTGTGAATTGGTGCTAAATAAATAAACTGAATGAAAAATAACTCAC

It looks like my current TRF output is sorted with the chromosome names at the top, and the microsatellites found at each chromosome following underneath the name. I think I can include the appropriate chromosome names into each row with no problem, and I can easily rearrange the data into the appropriate columns. My main issue is with the information itself. As you can see, my current output includes data I don't need: period size, size of consensus pattern, percent of matches between adjacent copies, percent of indels between adjacent copies, overall alignment score, and entropy are not required for my desired output (unless some of these terms are synonymous with information included in the desired output and I am simply not aware). Furthermore, my output needs to include the minimum flanking on the left and right, repeat purity, default slippage, and default stutter. Which TRF options do I need to use for this information to be included in my output?


text:
Help with setting DimPlot UMAP output into a 2x3 grid in Seurat
On my merged seurat object of 6 samples, when I use the split.by function in tandem with the Dimplot/UMAP plot, all six samples are displayed in series along a commonly labeled 'UMAP_1' x-axis in an order of (these are arbitrary for simplicity, here) Z, Y, X, C, B, A. Instead, I would like to specify these plots to output as a 2 row x 3 column arrangement ordered (1st Row) A, B, C and (2nd Row) X, Y, Z.
I'm assuming I could use something like: group_levels <- c("A", "B", "C", "X", "Y", "Z") to define the order I'd like, but I'm lost on how to get this to arrange the UMAP plots into a 2x3 grid I'd like.
Appreciate any help!

can you post the code or show with an example dataset (e.g pbmc_small in Seurat) ? I think it is a bit hard to follow in words what you need

I am sure you would prefer an answer with working code, so i encourage you to do the same here

how to get this to arrange the UMAP plots into a 2x3 grid

You can use the ncol parameter, which will set the number of columns in the "combined plot".
And for future reference, if you get the help page of the function you are using with ?DimPlot, you will see the .... This means that you can feed additional arguments to other functions that DimPlot() uses. Apparently,  CombinePlots() is used under the hood and has an ncol parameter.

Thank you! ncol solved my problem! Here's the code I ended up using to get the columns and samples ordered, not sure if it was the most straightforward though, but works:
object$orig.ident <- factor(x=object$orig.ident, levels = c("A", "B", "C", "X", "Y", "Z"))
DimPlot(object, reduction = "umap", split.by = "orig.ident", ncol =3)

source: Stack Exchange

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

sort: votes


text:
Seurat FindMarkers() output interpretation
I am using FindMarkers() between 2 groups of cells, my results are listed but i'm having hard time in choosing the right markers. Do I choose according to both the p-values or just one of them? If one of them is good enough, which one should I prefer? I am interested in the marker-genes that are differentiating the groups, so what are the parameters i should look for?
             p_val avg_logFC pct.1 pct.2  p_val_adj
UBD   3.080608e-06 1.7755753     1 0.000 0.07683959
CFB   1.262305e-05 1.3989233     1 0.067 0.31485675
RGS13 1.548593e-05 0.9009480     1 0.200 0.38626551
PYDC1 1.896309e-05 0.9622537     1 0.133 0.47299636
BIRC3 4.099414e-05 1.0129472     1 1.000 1.00000000
ICAM1 4.379895e-05 0.8219610     1 0.200 1.00000000

as you can see, p-value seems significant, however the adjusted p-value is not. and when i performed the test i got this warning In wilcox.test.default(x = c(BC03LN_05 = 0.249819542916203,  ... : cannot compute exact p-value with ties
I am completely new to this field, and more importantly to mathematics. Please help me understand in an easy way.

You need to look at adjusted p values only. I'm a little surprised that the difference is not significant when that gene is expressed in 100% vs 0%, but if everything is right, you should trust the math that the difference is not statically significant. Why do you have so few cells with so many reads? Is this really single cell data?

The p-values are not very very significant, so the adj. p-value. You need to plot the gene counts and see why it is the case. It could be because they are captured/expressed only in very very few cells. VlnPlot or FeaturePlot functions should help. You haven't shown the TSNE/UMAP plots of the two clusters, so its hard to comment more.

I am working with 25 cells only, is that why? so without the adj p-value significance, the results aren't conclusive?

That is the purpose of statistical tests right ? At least if you plot the boxplots and show that there is a "suggestive" difference between cell-types but did not reach adj p-value thresholds, it might be still OK depending on the reviewers. But with out adj. p-values being significant and without seeing the data, I would assume its just noise.

I've added the featureplot in here. Is that enough to convince the readers?

For me its convincing, just that you don't have statistical power. Did you use wilcox test ? May be you could try something that is based on linear regression ? though you have very few data points.

yes i used the wilcox test.. anything else i should look into?

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


text:
How to get transcriptome FASTA file for viruses for Kallisto pseudo-alignment?
I guess the title is self-explanatory. I'd like to find the reference transcriptome files for a few human viruses. My RNA-seq samples are from human tumor cells that were infected with viruses. The main part of my analysis is quantifying human gene expression from these samples. However, I also would like to find what proportion of my raw data is originated from virus RNA.
So, I'd like to learn how to find the reference transcriptome data for human viruses. Moreover, I appreciate it if you could give me some insight into your workflow for approaching this problem. Finally, if there is any specific command line, R, or python tool for this specific task, I'd like to know about it.

Interesting question, presume you're looking for things like HPV. Its a bit sparse on detail of course, so beyond a one word answer I'm not sure we can help. What about just using Blat?

@Michael HPV is a close example but my specific use case is for engineered viruses for cancer-killing. I mostly use VSV, VV, HSV, or Ad. What do you using Blat? Does it have a command-line interface? Can I give it raw fastq files? Does it give me the proportion of virus originated transcripts in my human transcript samples?

You might be interested in the viral sequence data from CATCH. There are various different viral subsets available.
You could also try the NCBI virus website, filtering on human hosts. There's a "Download" link on this page that allows you to download the full sequence set:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?VirusLineage_ss=Viruses,%20taxid:10239&HostLineage_ss=humans,%20taxid:9605&SeqType_s=Nucleotide
You've asked a very specific question, so I'm giving you a very specific answer. It seems that what you actually want to do is different from the self-explanatory question.
For the purpose of creating an index for Kallisto or Salmon, it should be okay to use the full dataset. Those programs are designed to handle eukaryotic transcriptomes as input, which typically include many gene copies and isoforms, so transcript / genome duplication shouldn't be much of an issue.
Without knowing more about your particular application, it's difficult to provide additional advice. From what you've written, i.e. "what proportion of my raw data is originated from virus RNA", a hit to any viral sequence would be sufficient. For example, if you want any more specific information about viral proportion / species, then you should be using a metagenomic tool like Kraken2, rather than Kallisto / Salmon.

Could you please explain this CATCH data more? Can I get cDNA of a specific virus in FASTA with annotated CDS?
I already knew about NCBI virus database, the problem is that it doesn't tell you which genome/transcriptome is the reference one for many viruses and just gives you many different isolated strains. The second problem is that the FASTA file doesn't have separate head and sequence sections for each CDS in the virus genome but it has one head with all the sequences of the virus. So, I guess I can't use it for the quantification of viral transcripts with Kallisto. Am I wrong?

Thanks @gringer this is a great answer. In fact it is a superb answer that I'll be following as well.

I found a possible answer to my own question by reading this biostar page and also finding an NCBI page for downloading the reference genome file for viruses.

You should go to this NCBI page and search for your target virus.
Then click on the virus name and you will be redirected to a page like this which gives you an Entrez records table. You should click on the direct link in the genome row in the table.
Click on Genome Assembly and Annotation report  on the new page. You will be redirected to a page like this with some information on reference files submitted to NCBI. The last column in this table is for links to the NCBI FTP site. There are two options to access RefSeq or GenBank FTPs. You can choose either of them. I chose RefSeq. Then, you will be redirected to a page like this with different fasta files for that specific virus.
the file with the name GCF_002815995.1_ASM281599v1_cds_from_genomic.fna.gz is the transcriptome file annotated with CDS info which can be used for Kallisto alignment.

Note: If you couldn't find your desired virus genome in the NCBI page mentioned above, you can proceed with the method for generating your desired transcriptome file from genome and annotation file using gffread package. The details on this method are provided in the answer to the Biostar question mentioned above.


text:
Transcribing DNA exons and append transcripts into a list
I am trying to find a way to read and transcribe a list of DNA sequences (list of lists) only when the for loop finds a start codon (triplet of the list items) and until it finds a stop codon, over and over, but the if function will only transcribe the first codon and won't go through the whole sequence.
I'm trying to do this with a mix of normal python and biopython functions cause I think it's simpler.

This is an example:
Start codon: ATG
Stop codon: TAG
gene_1= 'ACGGACTATTC'
gene_2= 'GGCCATGAGTAACGCATAGGGCCC
gene_3=GGGCCCATGACGTACTAGGGGCCCATGCATTCATAG
So, the first gene does not contain any transcript, whereas the second contains 1 and the third contains 2. I'm trying to get rid of everything outside these reading frames and append these transcripts into a list that should look something like this.
mRNAs = ['AGUAACGCA', 'ACGUAC', 'CAUUCA']
This is what I have so far:
genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    for i in range(0, len(ORF), 3):
        codon= ORF[i:i+3]
        if codon != 'ATG':
            next(codon)
            if codon == 'ATG':
                mRNA=codon.transcribe()
                if codon == 'TAG':
                    break
    return(mRNA)

mRNAs=[]
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)
        mRNAs.append(mRNA)
print(mRNAs)

But it is not really giving anything back, I wonder if the code it's too redundant and I really don't need to define a function here, do you know any better way to do this?
Thaaanks!!

There are a lot of issues with this code, some very basic Python mistakes such as forgetting quotes around ATG/TAG in your if statements,  and looping through your genes per base, instead of per codon as your code seems to suggest you wish to. Try adding print statements to your code to see at which positions you are not getting expected results, and edit your code accordingly.

Is this a homework? If yes, you should mention it's a homework. If not, are you sure you want to do gene prediction like this?

I agree, I forgot the quotes in the start-stop codons, however, the question goes more about the methodology used here, until now, I just managed to transcribe the start codon, but it wouldn't keep extracting the other codons.
@KamilSJaron This is not a homework, it's just an example of what I want I'm trying to do here which is extract fragments within the start/stop frame

I bet my shoes this problem must have a reliable solution. What about https://www.ncbi.nlm.nih.gov/orffinder/ for example?

Debugging 101: add print statements to see what you are actually processing. You are calling mRNA= transcription(codon). What is the value of codon? Isn't it a single base? Isn't your function expecting an array instead? have you even tried running this? Don't you get an error message?

I'm afraid none of your code makes much sense. You have written a function that expects an array, and then you are passing it a single letter
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)

Then, the function also doesn't make sense. What does next(codon) mean? How would your script know what a codon is or how to move to the next one? And your if codon == 'ATG': is within the if codon != 'ATG': block, so it will never be executed by definition.
I think what you are looking for is something like this, using the same basic approach as your script:
#!/bin/python3
from Bio.Seq import Seq
from Bio.Seq import transcribe

genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    foundStart = False
    foundEnd = False
    for i in range(0, len(ORF), 3):
        codon= "".join(ORF[i:i+3])
        if codon == 'ATG' and not foundStart:
            foundStart = True
        if foundStart and not foundEnd:
            cc=transcribe(codon)
            mRNA = mRNA + transcribe(codon)
        if codon == 'TAG':
            foundEnd = True

    return(mRNA)

mRNAs=[]
for gene in genelist:
    mRNA = transcription(gene)
    mRNAs.append(mRNA)
print(mRNAs)

Which returns:
$ foo.py
['', '', 'AUGACGUACUAG']

Of course, it only returns one mRNA since you only have one gene with an ORF in the 1st frame and your script was only looking at the forward strand and the first frame. In reality, ORFs can be in any of the 6 possible frames. You are also assuming that all ORFs will end with a TAG which is perhaps the case in the species you are working in? I don't remember if there is any genetic code that has TAG as the only stop.
If you do actually need to find ORFs in any frame, then I hope the script above should be enough to get you started.

@EdwinOrtegaArzola you're welcome, but please read the last paragraph. There are some serious issues with this approach and you shouldn't be using it for any real work. I strongly urge you to use an already existing tool for this instead of writing your own.

This is precisely the structure and functions I was looking for!! I forgot about the True and False statements and they are extremely useful here! so Thank you!
I agree on the ORFs issue, so, I will work later on the different reading frames, as for the stop codons I was really only using TAG as an example, I am actually including other stop condons in the code. Once, again thank you!

Just in case this isn't a homework problem, there are already plenty of programs available to do this. One program I use is getorf from the EMBOSS toolkit:
https://emboss.bioinformatics.nl/cgi-bin/emboss/getorf

Hello, this is not a homework problem but thank you for the link, I can use it to compare results!


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

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

samfile.close()

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

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

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


text:
Modifying protein_scale in ProteinAnalysis in SeqUtil in Biopython?
Is there a way to give weight arguments to each of the residues when using protein_scale in the ProteinAnalysis module of Bio.SeqUtils.ProtParam? I thought this would be done with the window and edge arguments, but from the documentation and the question How to allow "half windows" in biopython? I understand that the weights are symmetric and calculated automatically given a window size.
I'd like to understand what the idea behind this was, just to check that I understand this correctly. For example, if one calculates eg. hydrophobicity of a protein, why should it be a value per residue that depends on a particular window?
I was looking for a functionality like where the value given for each residue in the ParamData dictionary would be multiplied by a specified weight, like so:
```Y = ProteinAnalysis('MSTAEFAQLL')
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10)
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,edge[0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.])```

Of course I can code this myself, but was wondering if that functionality already existed.

There is not a way to control the edge weights as written, but making a modified version is straightforward. Namely make a class that inherits the class Bio.SeqUtils.ProtParam.ProteinAnalysis and overwrites the method relevant methods.
Having a gander at the files one sees the following:
def protein_scale(self, param_dict, window, edge_list: list):
    """[...]
    """
    # generate the weights
    #   _weight_list returns only one tail. If the list should be
    #   [0.4,0.7,1.0,0.7,0.4] what you actually get from _weights_list
    #   is [0.4,0.7]. The correct calculation is done in the loop.
    weights = self._weight_list(window, edge)
    scores = []
    ...

So the weights are generated by the 'private' method ._weight_list. Underscore means it is for internal use and for production code it is bad to assume that the name of the private method many not change, but Biopython has not changed in years so you'd be safe. And this is for a quick hack anyway...
class JankyProteinAnalysis(Bio.SeqUtils.ProtParam.ProteinAnalysis):
    overriding_weights = None

    def _weight_list(self, window, edge):
          if self.overriding_weights is not None: # dodgy override!
               print("Weight method hack done!")
               return overriding_weights
          else: # behave normally
               return super()._weight_list(window, edge)

So you'd do:
Y = JankyProteinAnalysis('MSTAEFAQLL')
Y.overriding_weights = [0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.]
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,float('nan'))

Note that the edge argument is no longer doing anything.

Thanks so much! Super helpful.


text:
What are the symbols of $*$ and '$\_$' in an unknown alignment format for HLA data from the IMGT dataset?
Has anyone worked with the IMGT HLA database/dataset before? IMGT-HLA git repo they have some convenient text files (eg. link to gene A .txt file) with the genomic data for the different alleles of each gene; what is the name of this type of alignment format? sample:
# file: A_gen.txt
# date: 2020-07-13
# version: IPD-IMGT/HLA 3.41.0
# origin: http://hla.alleles.org/wmda/A_gen.txt
# repository: https://raw.githubusercontent.com/ANHIG/IMGTHLA/Latest/alignments/A_gen.txt
# author: WHO, Steven G. E. Marsh (steven.marsh@ucl.ac.uk)

 gDNA              -300
                   |
 A*01:01:01:01     CAGGAGCAGA GGGGTCAGGG CGAAGTCCCA GGGCCCCAGG CGTGGCTCTC AGGGTCTCAG GCCCCGAAGG CGGTGTATGG ATTGGGGAGT CCCAGCCTTG
 A*01:01:01:02N    ********** ********** ********** ********** ********** ********** ********** *********- ---------- ----------
 A*01:01:01:03     ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
 A*01:01:01:04     ********** ********** ********** ********** ********** ********** ********** ********** ********** **********
...

There is the .msf format of the same dataset at, HLA gene A in .msf format, git repo, and there are 2 types of characters in the data [A,C,G,T,.], where the '.' denotes the gap (correct?). For this above, unknown to me, format there are 3 extra symbols, [A,C,G,T,.,$*$,|,_]. I assume that the symbol '|' (pipe character) is for the exon, then what do the symbol characters $'*'$, '$\_$', and '.' denote in this format? I assume that the '.' still denotes the gap, so what is the $'*'$ and '$\_$'? Does $'*'$ reference missing data in terms of the reference sequence? If so the .msf format lacks this information? and would it be reasonable to assume that '$\_$' is that the data at that base pair is equal to the reference sequence as sequence 1?
If my assumption is correct and this new dataset accounts for 'missing genomic data' in the alignment and the msf does not, then the .msf is penalizing with a gap data that is actually missing?
(the FASTA data for this gene, link to the fasta file, shows that the reads for allele are of a different length)


text:
a proper Design Matrix for several drug treatments with both control negative and control positive
I have a dataset of RNA-seq samples for testing different drugs on the presence of another drug. One of my samples is the normal cells with no drugs (control negative) and another is the cells with just the base drug (control positive). All of the other samples are cells with the base drug and another testing drug. So, my question is how to do design a proper Design Matrix for comparing these samples. I can separately compare with control negative and control positive. But, my intention is to somehow neutralize the effect of one of the controls in all the samples and compare to the other one. Is that possible? How should I design the matrix?

Literally one undrugged control?

Yes, Literally! Unfortunately I don't have biological replicates for any of these samples.

That's unfortunate.  I guess you can do some ratios in excel, but you can't do anything sophisticated without replicates.  And your conclusions will be really shaky.  Wet lab people are going to have to test your conclusions with better experiments.

I know that. But the whole experiment's purpose was to generate some hypotheses for further testing with replicates. So, I guess it can give me some insights!

Or it might give who a whole lot of false positives that waste time and money.  Everything  hinges on your one replicate being perfectly representative of what the control should look like.

I totally agree and discussed this issue a lot in the lab. But people still come to me with their no replicate datasets and I have to give them some results or insight! What can I do?

Give them a big 'told you so' 2 years down the line :-)

Are there any technical replicates? Could you share a table or list of the amount of samples and their conditions?

Condition is as I described. One sample is Mock, another sample is Drug1 and all the other samples are Drug1 + Drug2-7. Deug2-7 means 6 different drugs which are being combined with drug1 separately to see their simultanous effect. What do you mean by the amount of samples? Like number of Reads?


text:
Differential Gene Expression with Replicates for some of the samples
[this question has also been posted on Biostars; some additional clarification from there has been copied into this question]
I've been asked to analyse a set of samples in which their control sample has been done in duplicate (biological replicate) but the other samples have just one sample each.
The analysis purpose is just hypothesis generation for further experiment design. So, is it possible to have duplicate samples for one or two of the groups and use one sample for other groups? Could we estimate a constant dispersion from these replicates and apply it to all of the samples for differential gene expression analysis.
I'm using EdgeR/Limma for my differential gene expression analysis. When there are no replicates the eBayes function gives me an error. However, if I add a replicate to just one of the samples, it won't give any errors. Does this mean that it uses the calculated dispersion from the duplicate condition for all of the other one sample conditions? Or is it a bug that it doesn't give an error?
Should I consider this method of using duplicates for some samples instead of just using no replicates in such large experiments that are being done for hypothesis generation?
I have read the EdgeR recommendations (page 23-24) on what to do without (enough) replicates, but the explanation is not complete. Like it doesn't say how to find these housekeeping genes. About its first case of using a common fixed dispersion. I guess what I'm asking here is I'd like a more reliable measure of this common fixed dispersion. I say instead of using a theoretical dispersion why not try to estimate dispersion from the few duplicate samples and then use it for all the samples. I'm asking if someone has already done this thing and I'd like to know about their experience and how reloable could thing method be.

About the affordability is more expensive to do several inconclusive experiments than one well done. The term you are looking for are technical replicates it is a know and old method used (check the limma vignette).

I don't mean techincal replicates. I meant biological replicates all along. But just part of my samples have biological replicates, others ones are done in one sample each.

@gringer this question asks if people had the same experience and how they handeld that. What I'm describing is what I'm working on now. There is a set of sample which their control sample has been done in duplicate but the other samples have just one sample each.


text:
Bonferroni count = 1 in GO-Slim overrepresentation test (PANTHER) is it possible?
I encountered a problem while performing analysis using the PANTHER website. I run the statistical overrepresentation test for a genelist of 99 genes (GO-Slim BP, Binomial, Bonferroni) in February 2020. My result was Bonferroni count = 1 and a list of several enriched GO terms. Next we run the same analysis with the same settings in July 2020 (two days ago) and the result is dramatically different, the Bonferroni count is 1770 and no enriched terms are found.
I am confused - is the result Bonferroni count = 1 possible? Why these results are so different?
This is how the correction is described on PANTHER website:

The expression data analysis statistics now include a Bonferroni correction for multiple testing. The Bonferroni correction is important because we are performing many statistical tests (one for each pathway, or each ontology term) at the same time. This correction multiplies the single-test P-value by the number of independent tests to obtain an expected error rate.
For pathways, we now correct the reported P-values by multiplying by the number of associated pathways with two or more genes. Some proteins participate in multiple pathways, so the tests are not completely independent of each other and the Bonferroni correction is conservative. For ontology terms, the simple Bonferroni correction becomes extremely conservative because parent (more general) and child (more specific) terms are not independent at all: any gene or protein associated with a child term is also associated with the parent (and grandparent, etc.) terms as well.
To estimate the number of independent tests for an ontology, we count the number of classes with at least two genes in the reference list that are annotated directly to that class (i.e. not indirectly via an annotation to a more specific subclass).

GO-Slim BP, Binomial, Bonferroni

This seems easy you because I suspect have run too many tests at 5%.
If you run 10 tests at 5% rejection threshold, the Bonferroni correction threshold is 0.5% - you divide the number of tests by the rejection threshold.
One possibility is the rejection threshold is given in number of genems: in one instance you've run a pairwise test with several thousand genes so the rejection threshold has been set accordingly. The other test you ran I suspect was a multi-variate test, thus the shift in rejection threshold measured in number of genes.

The altnative is that Bonf = 1 is in fact P = 1, in other words there is no significant difference at all.  If P = 1 then this could be no genes show deviation hence 1770, because 1770/1770 =1
Looking at it again the above seems a better explanation.

I am not sure if understand what you mean. I run binomial test with Bonferroni correction in both cases using the same gene set. Once in April, next in July.
I was using the PANTHER website so I have no access to the code but the settings and the input were exactly the same:
nr of genes: 99
statistics: Binomial
correction: Bonferroni

I've put an alternative explanation above, without knowing the basis of the test its the best guess.

The p-value is not 1, there are several terms with p< 0.05 in the first case (with Bonferroni count = 1)
Analysis Type: PANTHER Overrepresentation Test (Released 20200407)
Annotation Version and Release Date: PANTHER version 15.0 Released 2020-02-14
Analyzed List: Dupl_N1_WT_0503_list.txt (Mus musculus)
Reference List: Mus musculus (all genes in database)
Test Type: BINOMIAL
Correction: BONFERRONI
Bonferroni count: 1
(P-value)
1.19E-02
1.19E-02
1.19E-02

I am sorry for the formatting

No you ran the Binomial, Bonferroni. The red alert is the Bonferroni correction, because this is classically when things fall-over.


text:
What is the best way to address the question of doublets and multiplets in a single cell RNA seq data set?
I have attached a histogram plot of the number of genes per cell in a single cell RNA seq data set of lung endothelial cells. I do not find a bimodal or multimodal distribution of the number of genes per cell within any particular condition. Doublets/multiplets should have an aberrantly high number of gene expression and if many, can lead to a bi or multimodal histogram plot of the number of genes detected per cell.  Can I conclude that I have an insignificant presence of doublets or multiplets in my data set from this plot and hence continue with downstream analyses?
I am reluctant to use automatic doublet-finding algorithms because they have a high possibility to remove cells in a transitional differentiating phase. I also think that the default threshold of 2500 genes per cell given by Seurat might sometimes lead to the loss of good cells depending on the cell type.

Apart from trajectory analyses, are there any other analyses on single cell RNA seq data that are sensitve to doublets/multiplets?
Thanks in advance for your help.

You should not expect a multi-modal distribution. See this previous discussion: https://bioinformatics.stackexchange.com/questions/3211/determining-doublets-in-single-cell-rna-seq

Indeed there's no satisfactory doublet/multiplet identification methods as of 2020.
The most reliable way to remove doublets is still manual removal of doublet subclusters guided by marker genes.
Not only trajectory analysis can be impacted by doublets, but clustering and embedding can also be affected for obvious reasons.


text:
What is a good rule of thumb for the threshold of noise versus signal for RPK in RNA seq?
I have RPK values (RNA seq) and I'm wondering what is a good rule of thumb for what is considered to be noise versus what is considered to be signal? I.e what should I choose as a threshold value for what is considered noise versus not noise? (just as a rough guideline)

Signal and noise are best determined at the count level, rather than after normalisation. I don't think there's an easy way to get a signal level from normalised data.
When I wanted to establish the noise level of counts from Illumina reads from lots of different samples, I looked at the total count distribution for genes that had counts in only one sample, and used the maximum value of that (or possibly an elbow of that distribution) for my signal threshold.


text:
Recommended workflow for RNA long read Kallisto-like TPM estimation?
On short reads, Kallisto/Salmon is a standard workflow for measuring RNA transcript counts in TPM. However, when I tried to Google a similar workflow for long reads, it's not clear there is a replacement.

https://github.com/BrooksLabUCSC/flair was excellent. Highly recommended. One could also ask for TPM.

While the selective-alignment algorithm that is part of salmon is designed for short reads, salmon also supports quantification via alignments in a SAM / BAM file (aligned by other tools).  You have to be somewhat careful about setting the proper command line flags, but nonetheless, the quantification should work fine.  You can, for example, find a pipeline that uses minimap2 for alignment followed by salmon for quantification of long reads here.

Stringtie now does long reads too. We've found flair to be quite difficult to install and to crash often. That said, stringtie isn't very good about showing you its working (e.g. I have no idea how to get it to give the assembly it uses to figure out all the isoforms). Its values for TPM/FPKM also seem a little weird, in that they're perfectly correlated with each other, and for transcripts they're perfectly correlated with coverage too.
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1910-1
https://github.com/gpertea/stringtie
TALON is another tool that does isoform assembly. I'm not sure whether it does quantification though.
https://github.com/mortazavilab/TALON#


text:
How to understand methylation level?
For Bisulfite Sequencing data, we can calculate methylation level for a specific locus,
$$
\mbox{methylation level} = \frac{\mbox{# methelyated reads}}{\mbox{# total reads}}
$$
Here I think, for this locus, it can only have one status, methelyated or not. Or to say, methylation level here can only be 0 or 1. Why does it have other values in (0,1)?

You are right that a single molecule in a single position is either methylated or not methylated. However:
First, assuming your organism of interest is diploid (or of higher ploidy) one of the chromosomes could be methylated, the other not. That would give you a level of 0.5 and can be found in imprinted regions (where the paternally inherited chromosome is methylated and maternal not, or vice versa).
Second, except if you use a pure cell population of a single cell type it is likely that there will be cellular heterogeneity in your sample, e.g. from different types of blood cells or neurons vs glial cells in the brain. Some of these may have methylated CpGs in this locus (depending on cell-type-specific gene activity) while others can be unmethylated. Depending on the relative proportions of these cell types the methylation level can be anything between 0 and 1 for a mixture of cell types.
Thirdly, sometimes CpGs get grouped together. For example, 5 nearby CpGs are combined to give a single methylation level for this locus, and then the level can be 0; 0.2; 0.4; 0.6; 0.8 and 1.
Lastly, there could be technical artifacts or technical reasons why a position was not called as methylated (insufficient quality for a single read?), leading to deviation from 0 or 1.

Thank you very much! I noticed some databases supply a single methylation level for each promoter region, which I could not completely understand before. Now I think that's the third case in your answer. Your answer is really helpful.

Thanks for the explaination. One further question though, why would the methylation distribution are biomodal when 1) methylation could be heterozygous and 2) we have mixture of cell type where each type represent different methylation pattern?


text:
Bootstrapping cell type identifications (scRNA-seq)
I have scRNA-seq data held in a SEURAT object, and cell-types were identified with SingleR (essentially a vector of strings - Monocytes, Granulocytes etc.). I have two conditions - treated and untreated. I want to know if the difference in the relative frequency (%) of a given cell type between the conditions is significant. For this, I thought to perform bootstrapping, to obtain frequencies for many different 'simulated' datasets - and then perform a statistical test. However, I cannot understand how best to do this. Should I subsample the scRNA-seq data, with replacement and calculate the frequencies each time? Can somebody provide a step by step workflow on how to do this?

I would use permutation:
i) Get the difference that you are interested in, i.e. number/proportion of monocytes in one group vs the other.
ii) Get the relevant columns (condition, cell label, ...) from the @meta.data slot of your seurat objects and rbind().
iii) Shuffle the cell labels randomly and split the data frame into two (per condition).
iv) Compare the difference between the cell counts/proportions of interest from the two permuted data frames and save this value.
iv) Repeat iii and iv for 100s of times and calculate how many times have you observed a difference at least as extreme as your real observation (i). Once divided by the number of trials, this would give you a p-value stating how likely it is that you get the observed value.
I have to add that batch effects are severe in scRNA-seq and hence a linear model where you might be able to account for batches could be more helpful.


text:
Using variables with fasterq-dump?
I am trying to download multiple fastq files from the SRA NCBI database which, conveniently enough, have their IDs in a range. Building a simple loop doesn't seem to do the trick, as bash appears to interpret my variable as an additional command? Help would be appreciated. This is what I have for now:
for i in {256..345};do
    sra_code = "DRR163""$i"
fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8

done
```

You have an obvious error in your shell script that should have given you an error message:
sra_code = "DRR163""$i"

Should be:
sra_code="DRR163""$i"

You cannot have spaces around the = sign in a variable assignment. The line sra_code = "DRR163""$i" means "run the command sra_code with the arguments = "DRR163""$i". Since there is no such command, you see an error.
Here's a working version of your command:
for i in {256..345}; do
    sra_code="DRR163$i"
    fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8
done

which, conveniently enough, have their IDs in a range

I think it's a bit risky to rely on that - for this particular case they are not going to change but you may need to download more samples in the future.
I would suggest avoiding sra and fastq(er)-dump altogether and use ena instead.
In your case, it seems you want the runs in project PRJDB7736. From here https://www.ebi.ac.uk/ena/browser/view/PRJDB7736 you can get a tabular file (see link in Download report ... TSV) with the sample characteristics. Column fastq_ftp gives you the ftp link to the fastq files that you can pass to wget or curl.
You can also get the same report programmatically with:
curl -s -X GET "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJDB7736&fields=run_accession,sample_description,library_name,library_layout,fastq_ftp&result=read_run" > ena.tsv

You should edit the fields=... list to suite your needs. The fields available are those from:
curl -X GET "https://www.ebi.ac.uk/ena/portal/api/returnFields?dataPortal=ena&result=read_run"

I don't know what you think... I find this route much more useful and transparent than sra's.


text:
Error running GenomicVis vcf.venn(vcf.files, 'GRCh38', sample.names)
I am trying to draw venn diagrams with a GenomicVis package (R).
library(GenomicVis)

It is loaded with a warning:
Warning messages:
1: replacing previous import ‘IRanges::shift’ by ‘data.table::shift’ when loading ‘GenomicVis’
2: replacing previous import ‘S4Vectors::second’ by ‘data.table::second’ when loading ‘GenomicVis’
3: replacing previous import ‘S4Vectors::first’ by ‘data.table::first’ when loading ‘GenomicVis’

f2 <- system.file("extdata", "S038_10_wes_S43.vcf", package ="GenomicVis")
f1 <- system.file("extdata", "S038_1_wes_S25.vcf.gz",package = "GenomicVis")
vcf.files <- c(f1, f2)
sample.names <- c('Sample 1','Sample 10')
v <- vcf.venn(vcf.files, 'GRCh38', sample.names)

But then this Error comes up, no clue how to solve it:
Error in rowData(read.vcf(f, genome, ...)) :
 could not find function "rowData

I also think that the problem is reading my vcf files as calling f1 or f2 results in [1] "", however when I open them with cat I can see vcf table as it is. The first line is ##fileformat=VCFv4.2. I tried using fullpath or running R from this folder and still the file was not loaded as system file is used only for loading example files from a package.
Loading file with a readVcf('S038_10_wes_043.vcf', 'hg38') worked well.

Hi @EkatSim this might be enough for a member to trouble-shoot, however this could be a result of your input format. You might need to supply more details. Personally, I'd like to peer into f1 and f2 objects to see what is there.


text:
Power calculation for GWAS/EWAS
I want to investigate, how much sample size i needed to obtain 80% power for GWAS/EWAS studies. Phynotype trait is discrete (not case/control) for human disease.
I wonder, does anyone has came across this? Really appreciate any suggestion or any tool/software which can answer this?

Might be worth having a look at this website: http://csg.sph.umich.edu/abecasis/gas_power_calculator/

thanks, but this if for case/control only. Our study is NOT case/control. Any other suggestion?

OK - could you edit your question to include some more details about the study you are planning (i.e. quantitative/discrete traits, species, intended outcome etc) and it might make it easier for people to give a good answer to your question. cheers.

I just Edited the question.

Do you mean continuous? If it's discrete and there are only two conditions, then the power calculations should be the same as for a case/control study.

@gringer yes I was also confused by this. what is the actual trait you are looking at?

Sorry for confusing, its a trait in young girls about their early stage in sexual maturation that happens in some higher primates and in humans, typically peaks at around 20 years of age, and is eventually involved in the development of pubic hair, body odor etc.

OK, so the trait is a continuous variable? In this case, it seems like it can be a 'time' value, like 5 years, 6 years, 5.1 years etc?

Interesting question I would recommend the effective population size $N_{e}$,  and error bars around this. Ne can be calculated via $F_{st}$ and is perhaps the most important parameter to determine the genetic population size.
Under the Fisher definition $N_{e}$ is the number of breeding individuals in the population.
The central issue with the approach you are suggesting is that population size does not equate to genetic population size. Thus the population could be huge, but if it is genetically homogeneous it will appear very small. In contrast the number of individuals could be very small, but the genetic variation between between them very large so $N_{e}$ would be huge.

Appreciated promt response.            Actually, our population  is humogeneous, and have very low sample size.  I wonder, is there any (genetics related) software where i can play around Ne and Fst?

Usual suspects Arlequin etc.. because you have low diversity the traditional stuff will work


text:
why in RNA seq don't we only use reference transcriptome?
I would like to ask why in RNA seq analysis (alignment step) we use sometimes reference genome instead of reference transcriptome? thank you!

Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.
But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.

thank you! Could I use STAR or a similar program (not Kallisto or Salmon etc) with reference transcriptome?

You can make a STAR reference out of a transcriptome, but why would you do that?  It wasn't designed to be run like that.

if i am only interested in differential expression, could i use star with reference transcriptome to decrease computational effort?

As a rule, do not use programs for things they were not designed for unless you are absolutely sure you understand what you are doing and why.  If you want something faster than STAR, use Kallisto.

I actually think you could run STAR on a transcriptome. In their manuscript, you can run alignment using a GTF file.

I didn't say you couldn't, just that there's no good reason to.  The point of a GTF file is to mark features, especially exons,  in a large genome, so there's little point in using one on a transcriptome.  And of course you can't align to a GTF alone.

I'm confused, the GTF file is the reference transcriptome no? And the FASTA file the reference genome (of course both are needed to map to transcriptome)

@Chris_Rands I wouldn't have spontaneously called a GTF file a transcriptome. For me, a transcriptome would rather be a set of sequences (the sequences of the potential transcripts). But maybe this is not the common usage.

@bli I agree, the transcriptome is maybe a set of mRNAs and not simply a tabular file format of course, but I was thinking in terms of how this information is stored

@Chris_Rands I never worked with a transcriptome, but I imagined that this was stored in the form of a fasta file containing all possible transcripts. And I imagined that mapping on a transcriptome was done via an index built on this fasta file. Is it actually more common to map on a genome and then infer transcripts based on information present in a GTF file?

@bli i think splice aware aligners need the full genome FASTA but others use transcripts FASTA


text:
Extracting all information about a sample when using xtract from e-utilities
I would like to extract all information about each SAMPLE after running the following query (run the query and add a | grep SAMPLE for clarification on what I mean by SAMPLE):
esearch -db sra -query PRJNA514750 | efetch -format xml
I tried the following:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern EXPERIMENT -element SAMPLE
but this returns nothing (PS: SAMPLEs are within an EXPERIMENT tag). I read in the e-utilities guide that -pattern will divide the data into rows, and -element into columns, so I'm presuming that this didn't work because SAMPLE has multiple tags within it. So I then tried:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element random_SAMPLE_tag

where random_SAMPLE_tag is any tag within SAMPLE.
Here's a concrete example:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element TITLE

This works, but I want to get all the information about each SAMPLE, and I do not know beforehand what the tags within it are (I manually got TITLE in this case), and since I want to get this info for a quite a few studies, I can't manually check this.

You can give pysradb a shot.
study_accession experiment_accession experiment_title                                                   experiment_desc
                                      organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument           total_spots total_size  run_accession run_total_spots run_total_bases
SRP178555       SRX5243234             GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq    GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245909                    Illumina HiSeq 2500  19052916     446793895  SRR8435699    19052916         971698716
SRP178555       SRX5243233             GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq    GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245908                    Illumina HiSeq 2500  18156224     428457264  SRR8435698    18156224         925967424
.. [truncated]


text:
Does the kinship and inbreeding coefficients depend on population frequency of an allele?
I am reading Section 5.2, Kinship and Inbreeding Coefficients, of Kenneth Lange, Mathematical and Statistical Methods for Genetic Analysis. There the kinship coefficient $\Phi_{i,j}$ is defined for two relatives $i$ and $j$ as the probability that a gene selected randomly from $i$ and a gene selected randomly from the same autosomal locus of $j$ are identical by descent.
I would think $\Phi_{i,j}$ should depend on the particular population frequency of the gene or allele at that locus. However, the book does not seem to indicate that or at least does not stress that at all. Is there a dependence or not? The same question applies to the inbreeding coefficient.
More importantly, I doubt this definition is mathematically rigorous.

Here is an example to make what puzzles me clearer. It must be that the probability I am computing is not what the kinship coefficient. But this is the definition means to me. Some expert please elucidate the correct meaning of the definition.
Consider the simplest pedigree tree of one family of two parents with brothers $A$ and $B$. Allele $a$ is detected at a particular locus. Let $b$ be the other allele of the same gene. We ask for the probability of the brother $B$ having allele $a$ conditioned on $A$ having allele $a$.
In the following, to clutter symbols, we abuse the symbols by use $a$ to denote the population frequency or the unconditional probability of finding allele $a$ at the locus for the whole population as $a$. The same goes for allele $b$. $a,b\in[0,1]$ and $a+b=1$. We solve this problem by listing as follows all possible parental gene configurations, the configuration probability (considering the symmetry of paternal and maternal loci), as well as the probability finding allele $a$ for either $A$ or $B$ at the locus.
\begin{align}
\text{parent config}~~~~ & \text{Pr(config)} & \text{Pr($a$|config)} \\
aa|aa~~~~~~~~~  &  ~~a^4   & 1 \\
aa|ab~~~~~~~~~  &  4a^3b   & 1 \\
aa|bb~~~~~~~~~  &  2a^2b^2 & 1 \\
ab|ab~~~~~~~~~  &  4a^2b^2 & \frac34 \\
ab|bb~~~~~~~~~  &  4ab^3   & \frac12 \\
bb|bb~~~~~~~~~  &  ~~~~b^4 & 0 \\
\end{align}
The joint conditional probability of both $A$ and $B$ have allele $a$ at the locus is Pr($a$|config)$^2$. So the required conditional probability is
$$P(r):=\text{Pr}(B\text{ has } a|A\text{ has }a)=\frac{\sum_\text{config} \text{Pr(config)Pr(a|config)}^2}{\sum_\text{config} \text{Pr(config)Pr(a|config)}} = \frac{r^3+4r^2+2r+4r(\frac34)^2+\frac4{2^2}}{r^3+4r^2+2r+4r\frac34+\frac42}, \quad r:=\frac ab\in [0,\infty).$$
$P(r)$ depends on $r$ which depends on the population frequency of the alleles as I have claimed.
Consider $[r^3+4r^2,2r,3r, 2]$ as a positive entry vector. $P$ is a convex combination of $v:=\big[1,1,1,\frac34,\frac12\big]$. So $P$ is between the minimum and maximum of the entries of $v$, and
$$\frac12\le P\le 1.$$
We see the extrema are achieved at
$$P(r)=\begin{cases}
\frac12, \quad r=0 \\[0.7ex]
1,\quad r\to\infty
\end{cases}$$

Interesting question. The question is what do you mean by a "population frequency"?

@Michael: I have added a detailed example to illustrate my question. Please review.

The definition as set by Kenneth Lange's book is indeed quite vague and thus unrigorous. It is a probability not conditioned on observing any particular allele at a locus thus not on the population frequency of that allele, but conditioned on the following where we focus on one particular locus. It is rigorously defined as follows
A graph of lineage contains nodes, which we call persons, of a particular locus. Each node, or person, contains exactly two genes. We draw one and only one directed edge (parental relationship) from a gene of a parent person to the gene of a child person if the latter is inherited from the former.  This is a directed acyclic graph (DAG). Each gene can be connected by no more than $1$ edge directed to (as opposed to away from) it. The two genes in one person cannot be connected by edges. Two persons are said to be connected if there is an edge pointing from a gene of one node to a gene of the other node. Two persons cannot be connected by more than one edge. For a pair of connected persons and their directed edge between them, the edge is equally probable to connect to one gene as the other gene of a person, i.e. the probability is $\frac12$. A pair of genes (persons) are said to be connected if there exists an undirected (simple) path from one gene (person) to the other. Two thus connected genes are called identical by descent.
In this setting, the DAG with respect to the persons is given and deterministic whilst any directed edge in the DAG assigned to two persons are random with independent probability $\frac12$ with respect to the genes within each person.
For a pair of persons on a lineage DAG that is connected, the probability of one gene equally probably chosen amongst the two genes of one person being connected, or identical by descent, to one gene equally probably chosen amongst the two genes of the other person is the coefficient of kinship. It can be proven that the probability is $\displaystyle\sum_{p\in P} \frac1{2^{n(p)}}$ where $P$ is the set of all distinct (simple) paths connecting one person (not gene) to to the other and $n(p)$ is the number of persons on person connecting path $p$.
As an additional assumption, the genes that are connected have to be the same and thus assume the same allele.

With regard to the conditional probability example in the question, the coefficient of kinship is not adequate to be used to compute the desired probability. Part of the reason is the question inquires about the existence of an allele rather than picking one gene with 1/2 of the probability. We need some more fundamental probability to compute that. We will demonstrate that as follows.

We will compute the probability of the two siblings having both allele $a$ as follows.

For a pair of full siblings, they can have exactly $0, 1, 2$ identical by descent genes. Let the respective probabilities be $r_0, r_1, r_2$. Notice we are not picking genes of a person at random but talking about the existence of particular types of genes in a person. For two full siblings, some thoughts lead to $r_0=\frac14, r_1=\frac12, r_2=\frac14$.
By symmetry and without the loss of generality, we fix the connecting genes in the two siblings of the respective edges. We only need to consider the connections of the edges relative to the genes of the parents.
(1) Conditioned on there is no identical by descent genes, the sibling genes can independently assume any allele. The probability of having the configuration of $ab|ab$ (not considering gene ordering) is $2^2a^2b^2$ considering the distinct ordering of the alleles, that of $aa|ab$ is $\binom41 a^3b$, and that of $aa|aa$ is $a^4$. The total probability is thus
$$Pr_0:=\frac14(4a^2b^2+4a^3b+a^4).$$
(2) When there is exactly $1$ pair of identical by descent genes, the two identical by descent genes can only assume the same allele, while the other genes can assume alleles independently. When the identical by descent genes assume allele $a$ the probability of the two siblings both having allele $a$ is $a$. When the two identical by descent genes assume allele $b$ the desired probability is then $ba^2$. The total probability is then
$$Pr_1:=\frac12(a+ba^2).$$
(3) When there are exactly $2$ pairs of identical by descent genes, each pair can only assume the same allele but independently of each other. The only way that the two siblings not to both have allele $a$ is for the two identical by descent pairs to assume allele $b$. The total desired probability for the two siblings to both have allele $a$ is then
$$Pr_2:=\frac14(1-b^2)$$

Now we compute the probability of any one of the siblings having allele $a$. Each gene can of course assume alleles independently. So the probability is simply
$$Pr_3:=1-b^2.$$

Therefor the desired probability of both siblings having allele $a$ given that one of them has allele $a$ is
$$Pr(a):=\frac{Pr_0+Pr_1+Pr_2}{Pr_3}=\frac{\frac14a(a^3-6a^2+5a+4)}{a(2-a)}$$
since $b=1-a$. Some tedious algebra which can be carried out by, say, Mathematica shows this is identical to the expression in the question in terms of $r=\frac{a}{1-a}$.
$$Pr'(a)=\frac12\Big(2-\frac1{(2-a)^2}-a\Big).$$
It is easy to see $2Pr'(a)-2$ strictly decreases for $a\in[0,2)$ as the last two terms decreases in that interval. $Pr'(1)=0$ and thus $Pr'(a)>0,\ \forall a\in[0,1)$. $Pr(a)$ thus strictly increases for $a\in[0,1]$ and $\frac12=Pr(0)\le Pr(a)\le Pr(1)=1$.

Thanks for your update

@Michael: You are welcome. The definition presented in Kenneth Lange's book is vague and unrigorous.

@Michael: My previous answer was still unclear. Now I have completely understood what identical by descent means along with the definition of the coefficient of kinship. I have revamped my answer. The description can be even more rigorous and detailed. But I want to keep it simple. Check out my answer and let me know if there is any point confusing.

Kinship selection has a long history in evolutionary genetics, it is where the whole "selfish gene" hypothesis was formed (Maynard Smith, Bill Hamilton, and George Price ). I would be amazed if it was not rigorous. The problem the inuitive understanding of how the co-efficient works in real data appears lacking. So for example $F_{st}$ simulation studies demonstrate that strong values occur above 0.2.
In terms of population size impact on the calculation, this is complex because it is a consequence of $N_{e}$, studies have been done on this - and in general extremely small populations don't adhere and tend to exhibit 'Muller's racket" mechanisms. The general opinion, e.g. in HW equilibrium, is that deviation from expected behaviour can occur due to a variety of factors not represented in the original equation, but given a reasonable  $N_{e}$, $F_{st}$ and $F_{is}$ the essential.
In every day terms the parameters can be explored by sensitivity analysis, so e.g. drive the genetic variation to extremes and look at the change in the coefficient, reduce the genetic variation and assess again. Look for past simulation studies etc ...

I do not see how this answers my question, though. I do not understand how this concept is defined under the ideal condition. Specifically, what is wrong with my example concerning siblings? Is my formulation of the definition of the kinship coefficient in the sibling case correct?


text:
why in 3' quantseq should we use reference genome and not transcriptome?
I am studying about RNA seq. I have found that in 3'prime rna seq (quantseq, tagseq) we should use a reference  genome  with    good    annotation  (plus known UTRs). My question is why?
Why could not we use reference transcriptome?
Thank you!

You asked a similar question before https://bioinformatics.stackexchange.com/questions/14093/why-in-rna-seq-dont-we-only-use-reference-transcriptome and were told to use the genome. If you so badly want to align to the transcriptome then do that but stop asking us for advice. Asking something when you are already convinced of an answer is a waste of time. Repeating the question when you don't like the answer is just plain rude.

It would be great if you can elaborate on where you "found" this information.
In the transcriptome annotation you have a defined number of transcript per genes, but there might be others that are present in your sample and not in the annotation. For example you have treated the cells with some extreme condition that changes the 3' processing machinery.
With 3' quantseq, you are left with the read fragment that is at the end of your transcript. If the ends of these transcripts are not part of the sequences in the annotated transcriptome, these will be unaligned.
These instances might be rare if your transcriptome annotation is good, but let's say if you find that a gene of interest which should be expressed is not, and you suspect it might be due to a different isoform used, you can always check it with the genome alignment (see where reads might be falling in the gene). With a transcriptome alignment, this will not be possible.

I'm not sure that's true:
https://www.nature.com/articles/s41598-019-55434-x
"QuantSeq. 3′ Sequencing combined with Salmon provides a fast, reliable approach for high throughput RNA expression analysis"

Yes I have seen this paper..thank you..but on another reliable source I have found the piece of information that I posted..


text:
Extract residue sequence from pdb file (in biopython, but open to recommendations)
I'm new to Biopython and I'd like to extract the sequence of residues from a pdb file.
My two questions are:

What is the simplest way to do this? (Esp. when there is more than one sequence)
and
Should I be worried about "PDBConstructionWarning: WARNING: Chain B is discontinuous"?

So far, I've obtained the residue sequence via:
p = PDBParser()
structure = p.get_structure("1ppi", "1ppi.pdb")
ppb=PPBuilder()
for pp in ppb.build_peptides(structure):
     print(pp.get_sequence())

seq = pp.get_sequence().__str__()

This seems to work well for this molecule. However, is there an easier way, especially when there is more than one sequence?
For example, I've read that one can also do
res_list = Bio.PDB.Selection.unfold_entities(structure, 'R')

but res_list is not a sequence of residues in str, and I don't know how to convert the output from res_list to a str sequence.
In addition (or perhaps because I'm going through the PPBuilder), I've recently gotten a lot of warnings of the type:
/usr/local/lib/python3.8/site-packages/Bio/PDB/StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line..
(For example, with structure = p.get_structure("5owu", "5owu.pdb"))
I've seen a discussion about how to silence such warnings, but should I be worried about this? I've also noticed that when I get these warnings, pp builder seems to give me more sequences then are there.

One thing to consider is whether you want the full sequence of the protein (equivalent to the SEQRES record) or the sequence of the residues with coordinate data in the file (from parsing the ATOM records). They often differ, especially for proteins with disordered regions, and the answers below cover both methods.

There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

Try help(PDBParser) for more, or see:
http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#__init__

The discontinuity is normal: it is a missing stretch due to lack of density. In PyMOL you see them as dotted lines.
I am sorry to say, the easiest/safest way to get the sequences is not using the PDB files...
import requests
data = requests.get(f'https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{code}').json()[code.lower()]
print(data[0]['sequence'])

One gets back a lot of information: example. This will tell you when you have multiple copies of the same peptide. The only information you do not get is what the PDB numbering is relative to the canonical Uniprot for that chain, which has to be taken from SIFTS.
A side note... Biopython.PDB is one of the best BioPython packages, but I personally I much prefer PyMOL as a python module (pymol2 module) to Biopython.PDB as it is functionally more complete and straightforward even if everything operates in a stage. You might want to check it out.
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1UBQ', 'prot')
    print(pymol.cmd.get_fastastr('prot'))

Thank you! I've accepted the other answer because of the point with how to iterate over the residues, but the info about the density/unitpro were very helpful, and I really appreciate the pointer to pymol!

You can also use the SeqIO module from Biopython https://biopython.org/docs/1.75/api/Bio.SeqIO.PdbIO.html
and do something like this:
from Bio import SeqIO

for record in SeqIO.parse(target, "pdb-atom"):
    print(record.seq)

where target is the path to your pdb file
the parser has multiple input/output functions, and pdb-atom is the one that will read out of the coordinates itself.

you can try PDBx Python Parser too:
http://mmcif.wwpdb.org/docs/sw-examples/python/html/index.html
it has an example showing exactly what you need
http://mmcif.wwpdb.org/docs/sw-examples/python/html/fasta.html


text:
DNA seq cutting using Linux
I need to use the cut command on a human insulin sequence row wise but i don't know the command.I already ran byte wise command but the gave me column wise results. I cant use the -f option since i need to use a delimiter with it and the seq does not have any delimiters.
please answer asap

Could you please show the commands you tried so far and how the result should looks like? If you could give a simple example of how your data, it would be a good help to understand your problem.

cut is a linux command, but might not do what you intend. Could you please describe your problem with greater clarity?

I have an assignment where i have to use Linux commands on a bioinformatics problem i choose  to cut sort and uniq the sequence. however the cut command worked for byte wise cutting but not field wise since fields dont exist in a seq. also can someone tell me why do we sort a seq because wouldn't that result in a mutation.

Use either generic *NIX cut (you do not need a delimiter if you use -c option) or head tools. Even better, use any of a number of bioinformatics tools for sequence manipulation, such as EMBOSS extractseq.
EXAMPLES:
# Write into file output.txt characters 1-20 from each line of file input.txt:
cut -c1-20 input.txt > output.txt

# Write the first 5 lines of input.txt:
head -n5 input.txt > output.txt

# It is easier if you have your sequence in fasta format.
# Write characters 1..20 and 100-110 from input.fasta to output.fasta:
extractseq input.fasta output.fasta -regions "1-20 100-110"

One of the easy ways to install EMBOSS suite is using conda (from bioconda channel).
SEE ALSO:
EMBOSS extractseq help: https://www.bioinformatics.nl/cgi-bin/emboss/help/extractseq


text:
Explaining the algorithm of RNA folding: what each symbol & value represents?
I need someone to explain to me from the nuts and bolts how algorithms/ maths is used to work out RNA folding. Explain it to me like I am an alien or child.
I am looking at this paper - https://eprint.ncl.ac.uk/240069
All the books on maths, and information theory & mathematical methods in bioinformatics are either too basic or fail to describe adequately what the research paper is doing in a way I can grasp.
I want to understand how ViennaRNA & NUPACK SPECIFICALLY uses maths to describe and work out RNA folding.
What do the math symbols mean, what are they representing, and why? Why do you write it like that?
In the paper I mention above, it describes how it calculates single stranded folding as per below?
a) What does all this mean?
x ∈ {0, 1}
|x |
= 1 if base i is specified in the
design of x to be bound and 0 otherwise. Further, let p
x ∈ [0, 1]|x |
^ I know conceptually it is referring to vectors and components, I know that the idea is to represent RNA folding as a graph structure in order to describe mathematically, but I do not know or understand exactly how or what it is describing in terms of the reality of an RNA folding. Two vertices connected by an edge represent a base-pair etc etc, that I get, but I cant seem to understand any further than that.
Here is the description from the paper:

Single stranded folding. Firstly, we evaluate the ability of a single
DNA strand to fold into its specified target structure (as shown in
Figure 1 top). Using the secondary structure predictor of the
ViennaRNA 2.0 so‰ware suite [13], we calculate the partition function
of all secondary structures that x might fold into.1 For any strand x,
let |x | denote its length and d x ∈ {0, 1} |x | a vector whose
component d x i
= 1 if base i is specified in the design of x to be bound and 0 otherwise. Further, let p x ∈ [0, 1]|x | denote a vector whose i-th
component p x i denotes the Boltzmann probability (obtained from the
partition function) that base i of strand x is paired with another
base. The single-stranded folding score Ssf is then defined as the
normalized Euclidean distance || · || between d x and p x as Ssf(x) =
1 − 1 |x | ||d x − p x ||. (1) Note that 0 ≤ Ssf ≤ 1 and Ssf(x) = 1 if
x folds unambiguously into its target structure.

Do you understand hybridisation and how that is relevant here?


text:
Installation of PRANK MSA in WLS Ubuntu 20.04 LTS
I want to install PRANK on the Windows 10 Linux subsystem (Ubuntu 20.04 LTS), I have followed the installation instructions to no avail.
wget http://wasabiapp.org/download/prank/prank.source.140603.tgz
tar xvzf prank.source.140603.tgz
cd prank-msa/src
make
./prank

[ cp prank ~/bin ]
[ sudo cp prank /usr/bin ]

I followed those instructions but got stuck with cd prank-msa/src I don't know why they want me to enter a directory that doesn't exist. I also followed this video on how to install linux files but it wasn't helpful. Any suggestions would be appreciated.

Their instructions are outdated and don't match what exists on their server. The problem you're running into is not at the cd step but even before that - at your first step even. Your wget fails because that file does not exist. Go to http://wasabiapp.org/download/ and see for yourself - there is no /prank/ directory in there.
I tried using curl, and it wrote the 404 page as a .tar.gz file, so if you did that, the tar step would error out.


text:
Within and between sample count normalization
recently I came across a situation where RNASeq sample quality was expresses trough DESeq sizeFactor. So authors reported quality of their samples with respect to some publicly available datasets as good/bad if the scaling factor of their samples was equal/higher/lower than 1. What is more, they computed it using FPKM's from Cufflinks and executed geometric mean calculation over genes (as the one done in DESeq after log transforming the raw counts). Underlying logic being that with FPKM the samples were normalized "within" and with the additional DESeq, normalization "between", thus ready ready for any downstream analysis (differential expression, AS, etc.)
I haven't seen something like that before nor the use of DESeq scaling factor for estimating the quality of a sample (all samples are from the same species/tissue/condition/etc. )
Can anyone comment on why is this ok/not ok because I cannot wrap my head around the logic (then again I am just trying to replicate the study, no value contribution on my end.)
Thank you

could you please link to the paper in question?

The approach you are describing seems very strange.
Crucially, the Vignette for DESeq2 states that the model only works correctly with unnormalized counts as input:

It is important to provide count matrices as input for DESeq2’s
statistical model (Love, Huber, and Anders 2014) to hold, as only the
count values allow assessing the measurement precision correctly. The
DESeq2 model internally corrects for library size, so transformed or
normalized values such as counts scaled by library size should not be
used as input.

Normalized counts like FPKM are not supported and should be converted to counts using the appropriate import funcion.
I don't really understand the use of the normalization factors for quality control because the number of reads in a sample does not correlate well with the quality. Maybe you could link to the paper in case there is something hidden in the experimental design that would explain this.
If you are interested in quality control a useful metric in my experience is the 3' sequencing bias that can be calculated using, for example, RNA-SeQC.

The only reasonable use of scaling factors in QC is to ensure that they're not too out there, since if they differ by around 10x or more you likely had a sample that failed (you would pick this up in a normal PCA though) and need to exclude a sample.

Agreed, but a simple look at the number of reads per sample would tell you the same.

It sometimes doesn't, since it's often the case that you still have a number of highly expressed genes but little else. But yes, often this will be obvious when looking at MultiQC output or similar.


text:
Blastp MSA to the same length
To generate aligned protein sequence truncate to the same length, I use blastp in biopython for MSA. I have multiple queries and one subject sequence. However, the alignment never returns the same length.
How can I align all sequences altogether, thus resulting in a set of sequences of the same length?

Please consider accepting/upvoting the answers below.

@Matteo Ferla wrote in the comments,

It sounds like you'd want a MSA of the whole set and not a set of pairwise ones. If global MSA is what you actually want, try ebi.ac.uk/Tools/msa/muscle. Then to trim the ends before phylogeny there are dozens of tools.

The reason is because Blast aligns locally, i.e. a window, while a global alignment is the whole entirety of the sequence

BLAST isn't appropriate for multiple sequence alignment.
BLAST generates pairwise alignments. Furthermore, the aligning regions shown in a BLAST report only include subsets of the hit sequences.
To create multiple sequence alignment of hits you want to:

run blast
extract the hit sequences
perform multiple-sequence-alignment (e.g., using MAFFT)

Extracting hit sequences can be done using the command-line. Check the docs for blastdbcmd.
If you prefer a point-and-click interface, you can also do this using my lab's SequenceServer BLAST software by simply selecting the hits that you want to retain, and in the Download section, clicking "FASTA of selected hits". Or if you want all of them, click "FASTA of all hits".

Please add a statement stating that your lab created this software that you're posting about on multiple threads.


text:
Issues with AutoDock Vina
I am trying to use AutoDock Vina to do docking however I am getting this error. I would really appreciate any help as I have not been able to resolve this error despite trying numerous times.

The error is a bash error from your wrapper script. Which is failing to find vina. What is the detail of line 7 and does it match where you'd expect to find vina (i.e. where you installed it)?

Also, I am not a Windows person, but I believe cygwin is highly depracated in favour of the windows 10 dev mode ubuntu kernel. You might want to check for newer options.

@+matteo ferla Here is the wrapper script. My co-worker was able to run the script but I performed the exact same steps as him but I couldn't run the script. Do you think AutoDock Vina should be installed or something. The issue seems to be like you stated the script is not able to detect vina.
Wrapper Script:
"#! /bin/bash
for f in ligand_*.pdbqt; do
b=basename $f .pdbqt
echo Processing ligand $b
mkdir -p $b
vina --config conf.txt --ligand $f --out ${b}/out.pdbqt --log ${b}/log.txt
done"

So you can either use the alias command to tell terminal what vina means, or you can add the vina binary to /usr/bin/.

@MatteoFerla will give that a try and get back to you. thank you.

@MatteoFerla i added the vina application files to the "bin" folder but I am getting the same error. Which file in the screenshot below is the vina binary you are referring to. Any additional recommendations on what might be causing the error.
Image: https://ibb.co/ZzMV2Vp

You can edit your original post to add extra data btw —including SE image upload. If you run the command ./vina in the bin folder it should work and it's just cygwin oddity —which best left alone.
If so, modify your script to call /usr/bin/vina as opposed to simply vina.

@MatteoFerla
I added the vina files and wrapper script to the usr/bin folder and retried on Linux but I am getting the same error. I added the screenshots to my original post above. Any suggestions.

.exe? That is a Windows binary, whereas you are trying in Linux, please download the appropriate one or change your bash script to a batch script and run in normal DOS terminal —for 6 ligands you don't need a for-loop though. However, however... may I suggest before embarking on a computational biochemistry task that you read around or watch a few videos about using the terminal?

@Matteo Ferla thank you, which videos would you suggest going through before embarking on the task.


text:
Time for running ADMIXTURE analysis
I am trying to run the ADMIXTURE software for the first time to analyse the structure of an in-house dataset of samples comparing to 1000 genomes project ancestral populations. I am trying to determine ancestry just from the X chromosome of those individuals.
To do that, I have downloaded the X chromosome file from 1000 genomes (ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.vcf.gz). I have removed duplicated with bcftools and converted to plink binary with the following command in plink 1.9:
plink --bcf ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.bcf --keep-allele-order --vcf-idspace-to _ --const-fid --allow-extra-chr 0 --split-x b37 no-fail --make-bed

Next, I modified the .bim files to change the polymorphisms ID so them where the same as their position, as a way to standarize their names. Then I filtered related individuals, also with plink and found the intersection of .bim files from 1000 genomes and my in-house dataset (a bed file with microarray genotyped variants only from X non PAR region), that was already filtered to contain only variants genotyped against the plus strand (without indels).
Then I merged both datasets with the following plink command:
plink --bfile filename1 --bmerge filename2.bed filename2.bim filename2.fam --make-bed --out filemerged

That command generated a file with all the SNPS that would have more than 3 alleles in the merged file so I removed all them from both datasets and merged them again. On that merged file I pruned SNPs by linkage disequilibrium:
plink --bfile cromXmerged2 --indep-pairwise 50 10 0.1

, assign sample sex on the .fam file and filtered individual with missingnes higher than 10%:
plink --file mergedfile --mind 0.1 --make-bed --out

I also set all heterozigous calls as missing with --set-hh-missing plink command.
With the final input I runned a supervised ADMIXTURE as following:
./admixture finalinput.bed 3 --haploid="male:23" --supervised -j4

Before running of course I created the finalinput.pop file in open office calc, by creating a list with the assigned population for each sample, with my in-house dataset left blanc. I just used african, european and east asian 1000 genomes samples because my own dataset has brazilian admixed samples. The final dataset has approximately 1500 SNPs and about 1600 individuals. My computer set up is CPU: i3-2370M with two cores, 4 threads@2.4 GHz, RAM 8 GB@1333 MHz and an SSD 120 GB Sata 3.0. OS Ubuntu 20.04 64 Bits.
The thing is that ADMIXTURE have been running for three days without showing any error or producing result yet. I would like to know if it is normal to this analysis to take so long or if it could be any problem with my dataset.
I would appreciate any suggestion or advice!!!
Thanks in advance,
Iriel

This should take a few minutes, not 3 days! Definitely something seems to have gone wrong here. Are you running it locally or on a remote server? ADMIXTURE produces an output as its running, so you should check that it's doing that

I am running it locally. And it did not produced any output at the moment. I tried a unsupervised analysis before with K = 3 to 10 and it took a few days to produce Q and P outputs for the first 2 K values. But I canceled that run because supervised analysis would be better for my goals.

Did it literally not produce any output at all? I have just had a quick check and run it, and it produces this almost immediately:
Parallel execution requested.  Will use 8 threads. Random seed: 43 Point estimation method: Block relaxation algorithm

Nothing, it keeps showing these lines:743116 (QN/Block)  Elapsed: 0.354 Loglikelihood: -1.1648e+06 (delta): -103.745
743117 (QN/Block)  Elapsed: 0.331 Loglikelihood: -1.16479e+06 (delta): 6.94675
743118 (QN/Block)  Elapsed: 0.327 Loglikelihood: -1.16472e+06 (delta): 67.0046

I mean, It did not produce any file yet.

Seems like it's having trouble converging. There is a post here (https://www.biostars.org/p/368113/) which suggests that using  haploid mode and multithreading might cause issues converging and that you should try running it without multithreading. Actually, your dataset is so small that I doubt you need to multithread anyway :)

Thanks for the reply! I will try without multithreading then.

no problem. I am interested to see if this solves the problem.

Yes! I took just 6 seconds to run!

Great. Should add the solution as a response to this and accept the answer so other people who had a similar problem can find the solution.

Apparently there is a bug in ADMIXTURE that prevents for converging when using at the same time the haploid mode and multithreading (flags --haploid and -j respectively). The problem is solved by avoiding using one of them.


text:
How to modify DNA evolution model to fit actual data?
I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).
I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms.
I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.
The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?

If mutation is a random process, why would you expect to reproduce the same result at the end of a stochastic simulation? Most simulation methods will not give you the exact same results as the training data; the "truth" is only one instantiation of many possibilities. Can you say anything about the methods you are using?

@MaximilianPress I agree, there will not be a 100% fit. Even 90% similarity is good enough.
Using a phylogenetic tree, every time there is a branch, I copy the genome and evolve it independently based on branch length. i.e. if a random number is less than branch length, the nucleotide changes (for each position)

What is the purpose of the sequence simulation?

@AnthonyGuterres see my new answer, it was going to be too hard to fit into a comment. Again, hard to know how to modify your model without knowing anything about your model.

Mutation is not a random process because there are many mutational biases, e.g. the transition vs. transversion ratio. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

INDELible http://abacus.gene.ucl.ac.uk/software/indelible/
Evolver https://www.drive5.com/evolver/ (not to be confused with the PAML package of the same name)

I decided that this was too involved for a comment, so I'm making it an answer.
It is still quite difficult to answer the question given the absence of methods explanation; but I think that there's some misunderstanding of phylogenetic modeling going on here, and I can at least address that. I'm going to try to break it down into its components, starting with the basic statistical idea of simulation and then going back to phylogenies.
Statistical reasoning
Let's say that I flip 10 coins (2 outcomes: head (H) and tails (T)), and I observe the sequence HTHTTTHHTT. I then fit a maximum likelihood model, and find that the best model for the data has
$Pr(H) = 0.4$,
$Pr(T) = 0.6$,
and that these two outcomes are mutually exhaustive:
$Pr(H) + Pr(T) = 1.0$
I can then use a binomial model to re-simulate data based on these parameters. I will observe various sequences, for example:
HTTTTTTHHT
TTTTHHTHHH
HTHHTHHTHH...
However, based purely on the parameters of the model, I don't expect every new simulation to recapitulate the exact same HTHTTTHHTT sequence, because I have reduced its information to some rate parameters for the simulation.
Each of those new sequences will yield model parameters ($Pr(H), Pr(T)$) that are similar to the training data, but they will not necessarily be particularly similar to the training data in sequence.
Phylogenetic modeling
Phylogenetic modeling works in the same way, such that you reduce the observations at the tips of the tree into some model with some parameters. Most of the time you are estimating e.g. a rate matrix of parameters for transition between states for a set of residues (in simplest case for DNA 4 nucleotide states and for proteins, 20 amino acid states usually). If you are working with gene sequences you are probably using the 4 nucleotide states or some derivation thereof.
Simulations usually work by starting with some fixed ancestral state at the root (which can take various states), and then evolving forward from there along each branch of the tree using the rate matrix in a continuous-time Markov chain to "replay the tape of life", by assuming that the model you've fit somewhat reflects truth.
While opinions differ somewhat on what to expect from "replaying the tape of life" (see that linked review for more), basically no one expects to get the same results as the input data, or even 90% similar to the input data. That would indicate that evolution is highly deterministic, which just isn't reliably true. Of course there are situations where such determinism happens to be the case, but they are relatively rare, and I don't think that sequence evolution is considered to be very deterministic at all.
It's a little different in that the tree defines a correlation structure of the tips, which constrains the simulation results to be somewhat more similar to the training data than we expect from coin-flipping. But there's still no reason to expect particularly similar results on the tips.
But you can also model discrete characters on a tree, which can be even simpler. For example, I used the R APE package to simulate the following character states (A or B for each tip) on a simulated tree of 10 tips with the following commands:
tr = rtree(10)
c_tr = rTraitDisc(tr)

I can then build a model for this character and this tree using the following:
rates = ace(c_tr, tr, type="discrete")

and use it to simulate new sets of states:
# A (index 1) is estimated to be ancestral state with 99% confidence
new_sim = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim, adj = c(1,0))

Look! Not a single observation of state B! Let's try it again:
new_sim2 = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim2, adj = c(1,0))

Ok, this is more similar, but there's still another state B on the tips.
Conclusion
Evolution is not very deterministic. The correlation structure of the tree makes the outcomes less random than independent trials as in coin-flipping, but it is not at all expected that simulation results from phylogenetic models will match the training observations well. The best that can be said is that each simulated dataset can hopefully be reduced to the same model parameters as the input data, but that's a very different statement from yielding the exact same outcomes.
If what you are looking for from simulation is to yield the exact same sequence states at the tips, then you are unlikely to get that from phylogenetic simulation. You might occasionally get pretty close, but that is not at all the expectation.


text:
Gap-fill assembly with PacBio reads
I would like to gap-fill (or correct) a mammalian-sized assembly for which we have BACs that have been sequenced with PacBio CLR data. I have seen there are pbm22/gcpp tools available for correcting assemblies, but looking at how they've been used, they seem to correct existing assemblies built with FALCON or similar.
Can pbmm2 and gcpp be used with a reference constructed without pacbio data? E.g. given a mammalian-sized reference and a .bam file of PacBio reads for a BAC, what's the best way to gap-fill the current reference? Thx

My understanding is that PBJelly is the tool intended to solve this problem.
This person seems to have used it starting with Illumina scaffolds. I have not personally used it but have heard that it works pretty well.
Edit: note that the PacBio ecosystem is not super easy to use from scratch. However, they have made pretty good efforts to facilitate setup using conda, as for instance documented here. This person shows a conda recipe for pbjelly.


text:
Hybrid approach in metabolomics (NMR and MS)
Is it currently common in metabolomics to combine two quantification methods, namely NMR (nuclear magnetic resonance) and MS (mass spectrometry, either liquid or gas)? If so, what would be profits from such an approach? Is it a promising direction or people are going to stay with a simple approach and they will be sticking to either of these two?

Hi and welcome, perhaps you could explain whats yuor investigation?

Hi Michael, I am just curious if running an NMR experiment will boost the insights that I may obtain from pure MS experiment. I have the feasibility to run both of them, I am just wondering if it makes sense.

It strikes me as more of a wet-lab question than bioinformatics, but its just my opinion

Not necessarily, wet-lab question would be, for instance, about applicability of a given method to the samples that I have. But I know the pros and cons of both methods. What is of my concern is the next step - whether the analysis utilizing both of these brings something more to the table in comparison with sole MS / NMR.