text:
How can I predict if the binding site of a DNA-binding protein is located on the promoter sequence of a gene?
Is there any computational tool or any other method to predict if the binding site of DNA binding protein is actually located on a promoter sequence?

Possible duplicate of How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you predicted binding sites: genomic intervals where those TFs bind, with some statistical significance.
https://bioinformatics.stackexchange.com/a/2491/776
You can then map or intersect promoter regions with those TF binding sites with BEDOPS bedmap or bedops, respectively:
https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html
You can define your promoter regions however you like, e.g., from the TSS to 2.5kb upstream of the TSS. It's up to you and your experiment.
The use of bedmap is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.


text:
Quality GL entries in vcf file
I have a vcf v4.2 file with structural variants, the 8 columns describing the variant are clear, but in the genotype information I got confused by values of GL entries. In the vcf specification they mention that GL are likelihoods of the individual genotypes, but as far as I am concerned, likelihood should be a number between 0 and 1, therefore I am suprised to see negative numbers. Then I thought that it might be log-likelihood, but then it's strange that the higher the number is the better. Here is a sample of 5 variants with various called genotypes
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-22.4949,0,-4.49485:45:PASS:20:1:11:0:1:4:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-3.1888,0,-43.9892:32:PASS:16:26:13:2:9:2:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-0.190635,0,-49.2888:4:LowQual:26:51:31:2:10:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/0:0,-1.01435,-75.899:10:LowQual:18:54:26:2:14:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   1/1:-19.198,-1.5031,0:15:PASS:70:0:58:0:0:5:0:0

The vcf file was generated by delly genotyper. I followed their recommendations for germline SV calling. However, negative likelihoods are even in the vcf specification, therefore it seems that it's not something delly specific.
So, what exactly is GL in vcf files? And how do they relate to other entries GP (genotype posterior probabilities, not in the file by delly) and GQ (conditional genotype quality)?

The VCF spec clearly states that these are indeed log likelihood values:

• GL : genotype likelihoods comprised of comma separated floating point log10-scaled likelihoods for all possible
  genotypes given the set of alleles defined in the REF and ALT fields. In presence of the GT field the same
  ploidy is expected and the canonical order is used; without GT field, diploidy is assumed. If A is the allele in
  REF and B,C,... are the alleles as ordered in ALT, the ordering of genotypes for the likelihoods is given by:
  F(j/k) = (k*(k+1)/2)+j. In other words, for biallelic sites the ordering is: AA,AB,BB; for triallelic sites the
  ordering is: AA,AB,BB,AC,BC,CC, etc. For example: GT:GL 0/1:-323.03,-99.29,-802.53 (Floats)

With log-likelihood values, the higher the number, the better the value. So -2 is better than -3.

Dang, sometimes I get stuck in so complicated thoughts, while just the basic premise is wrong!

@KamilSJaron join the club! That happens to all of us :)

source: Stack Exchange

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

sort: votes


text:
how can i download rna-seq data sample for starting rna-seq analysis by R?
I am looking for an induction into RNA-seq analysis using R. I have yet to decide a data set, or a  subject area. Any help in this regard would be welcome.

Hi @Aminafzalifar welcome to the forum! My advice is have a think about what to focus on and return with a better defined question. Can you code, do you understand dataframes? Sauret and Bioconductor are popular packages in R which deal with this, however in my opinion (not universally shared) you need a local supervisor to help you.

I'm sorry, but we don't do this sort of open-ended discussion here. This site isn't a forum, it's strictly a questions & answers site and we only cover specific questions that can have specific answers. Your question is just too broad for this site's format. You might want to take the [tour] or visit our [help] to understand how the site works. Especially this page that explains the kind of question to avoid.

thanks and sorry about my bad question.

You must have found the airway dataset by googling "RNASeq R tutorial"
https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html
Why won't this work for you?

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


text:
Transcription factor targets for Human and Mouse genes
I want to get transcription factor targets for my list of genes (For Human and Mouse)
I have found ENCODE Transcription Factor Targets database (The link is given below). But I am looking for some other databases also if available. I will really appreciate if anyone can suggest any other database.
http://amp.pharm.mssm.edu/Harmonizome/dataset/ENCODE+Transcription+Factor+Targets
Thank you

What is the problem with ENCODE? I am aware of the existence of TRANSFAC, but I think is a non-free resource

@Plat My only problem with ENCODE is that I am not sure when was it last updated (I think in 2012)

JASPAR is another useful transcription factor profiling database.

The JASPAR CORE contains a curated, non-redundant set of profiles, derived from published and experimentally defined transcription factor binding sites for eukaryotes. It should be used, when seeking models for specific factors or structural classes, or if experimental evidence is paramount.

http://jaspar.genereg.net/

JASPAR is only based on TF motifs, which are very short sequences and not enough to precisely tell you the true TF-motif relationships, especially when multiple TFs have complex logic gates regulating overlapped targets.
An alternative way is to find some transcriptomic data upon perturbation of each TF. The drawbacks are that the redundancy of TFs will mask a lot of targets, and that the up/down regulated genes may not be direct targets but targets of targets (of targets of targets...).
Another way is to find ChIP-seq data for TFs you care about, and hopefully you can find the right biological context. But still ChIP signal isn't exactly binding pattern, which in turn isn't exactly indicating target relationship.
Overall, it's a challenging task to know what are the true targets of a TF.


text:
Looking for a good tool for repeat masking in long sequences
The repeat masker site has a size limit of 90kb.
Where can I mask repeats in a 500kb sequence?
Ideally, I need a thorough repeat masking of known and de-novo found repeats of various complexity. I am working on Windows but can access Unix too.

RepeatMasker can downloaded an run locally. There's then no size limit, though do note that it's quite slow for mammalian-sized sequences (you'll want to split it into chunks and run in parallel in such cases). For a small 500kb sequence that won't pose an issue.

RepeatFinder may suffice by Volfovsky et al. (2001) and is available for download.

From: http://meme-suite.org/doc/glam2_tut.html, here's a few tools that do various repeat masking tasks.

These tools mask various kinds of repeats: Repeatmasker, Censor, seg/pseg/nseg, DUST,
      XNU, DSR, GBA.

Censor, seg/pseg/nseg, dust, DSR all have locally available applications. DSR, GBA, and CAST are algorithms for repeat masking.

Please use this python program to mask elements -
dna = input("Please enter the DNA sequence or chromosome whose repetitive elements need to be masked:")
repele = ["aaaa", "tttt", "cccc", "gggg"]
counter=0
str1=str(dna)
for i in repele:
    str1 = str1.replace(repele[counter], "x")
    counter=counter+1
print(str1)


text:
install bioconductor
I used
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install()

for installing "bioconductor" but I encountered this error:
Error: invalid version specification ‘Bioconductor version cannot be validated; no internet connection?’
In addition: Warning messages:
1: In file(con, "r") : InternetOpenUrl failed: 'à¿q'
2: In file(con, "r") : InternetOpenUrl failed: 'pen the ð¥Xn'

Can anybody help me?

Hi and welcome to the site! Please [edit] your question and tell us: i) do you actually have a working internet connection? ii) What operating system are you using? iii) What version of R?

Hi. Thanks. 1- yes i have 2-  windows 7 3- 3.6.0 version

@Amina what is your local language and your keyboard configuration?. It might be a problem related to it. The command with this error is the install? Have you asked this in support.bioconductor.org (which is the official channel)?

@llrs I checked my system and everything was correct. I didn’t ask support.bioconductor. However, thank you for your advise.

I don't know what is correct for you, but definitely something is not correct, so please post your language configuration if you want further help. I suspect that part of the problem is from reading a string in a different encoding that the supported or something similar.


text:
Run Nextflow with file dependencies inside a Docker container?
I’m trying to run a Nextflow workflow in a custom Docker container. Without Docker, the workflow succeeds. But running it inside the container leads to an error because a dependent file cannot be found.
Here’s a minimal script to illustrate the issue (similar to the basic Nextflow script example):
params.reference = "/data/reference/GRCh38.fa"

reference = file(params.reference)
fastq_channel = Channel.fromFilePairs('/data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "$reference" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t "${reference}.bwt"
    """
}

Running this inside a container will cause the path of the query, /some/fastq/, to be mounted inside the container. However, the path to the database, /data/reference, will not be mounted. Consequently, I get the following execution error:

stat: cannot stat 'GRCh38-chr21.fa.bwt': No such file or directory

Without the container, the command succeeds.

I understand why the error occurs: Nextflow doesn’t track the dependency on $reference inside the command, and thus doesn’t mount its path in the container.
My question is, how do I fix this? The documentation is silent on this regard, and none of the workflows I can find online address the problem. For example, the GATK4 IARC workflows would presumably have the same issue, and aren’t configured for use with Docker.
Just adding the reference as a process input is insufficient, since BWA requires various other index files. Am I supposed to add these all manually as inputs?

The following workaround seems to work but ideally I don’t want to have to change the script portion of my process so that it works with Docker.
params.ref = 'data/reference/GRCh38.fa'

ref = file(params.ref)
ref_index_dir = ref.parent
fastq_channel = Channel.fromFilePairs('data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel
    file ref
    file ref_index_dir

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "reference/$ref" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t reference/${ref}.bwt && touch ${sample}.sam
    """
}

Here, reference inside the script section refers to the parent directory name of the reference, and needs to be adjusted accordingly. Having to refer to reference/$ref is annoying but at least it works — both inside a Docker container and outside of it.

Process input files need to be declared properly to allow nextflow to create the required docker mounts accordingly. Therefore the reference file should be declared in the process input block, eh:
input:
 set sample, file(fastq) from fastq_channel
 file reference

The issue is that BWA for instance doesn’t just use the “reference” file, otherwise there’d be no issue. But it additionally derives various other files from the reference filename opaquely (i.e. these names are in principle implementation details and not known, except that they all have the reference filename as a prefix).


text:
How do i map subcellular localizations to lower the number of classes to predict
For a project, I want to predict the subcellular localization of some proteins. I am a student in software Informatics at the University of applied Science HEIG-VD.
Subcellular localization prediction is a common problem in the bioinformatics, but I have a problem
With the data I found. In this database: of human and mouse proteins, we can see that there are more than 30 different localizations.
At my level, and for the use of my exercise, 30 classes to predict is too much.
Where can I find a mapping to reduce those 30 classes to 4 or 5?
I tried QuickGO, to see the ancestor tree of the different classes, but nothing conclusive. Somebody might have done that, but I lack knowledge on bioinformatic web research.
There are the classes I wanted to map to reduce the number of classes :

Cytoplasm (GO:0005737)
Cytoplasmic Vesicles (GO:0016023)
Endoplasmic Reticulum (GO:0005783)
Endosomes (GO:0005768)
Extracellular (GO:0005576)
Golgi Apparatus (GO:0005794)
Lipid Particles (GO:0005811)
Lysosomes (GO:0005764)
Melanosome (GO:0042470)
Mitochondria (GO:0005739)
Nucleus (GO:0005634)
Peroxisome (GO:0005777)
Plasma Membrane (GO:0005886)
Synaptic Vesicles (GO:0008021)
Cellular Component Unknown (GO:0008372)
Apical Plasma Membrane (GO:0016324)
Basolateral Plasma Membrane (GO:0016323)
Centrosome (GO:0005813)
Golgi Cis Cisterna (GO:0000137)
Cytoskeleton (GO:0005856)
Early Endosomes (GO:0005769)
ERGIC (GO:0005793)
Inner Mitochondrial Membrane (GO:0005743)
Late Endosomes (GO:0005770)
Medial-Golgi (GO:0005797)
Outer Mitochondrial Membrane (GO:0005741)
Secretory Granule (GO:0030141)
Golgi Trans Cisterna (GO:0000137)
Golgi Trans Face (GO:0005802)
Tight Junction (GO:0005923)
Transport Vesicle (GO:0030133)

One way to do this would be to plot the first 10 GO categories of the the CC sub ontology (CC from cellular component). Or you can just visit the official webpage to see the list of children terms of them.
However there aren't only 4 or 5 categories, so you could arbitrary pick some of those.

This is probably one of those things where as long as you can justify your choice somewhat.
It might be worth seeing what others have done (see loctree3 and PSORT) and seeing how well they match your case, then add in any subdomains you have a particular interest in.


text:
Extracting the number SNP in each range
I have called copy number and I have 2 files (I have shared the link of my files ); One contains some ranges
> head(cndata[,c(2,3,4)])
     Chromosome   Start      End
4470          1   51479   817980
4471          1  818499  1136753
4472          1 1138735  2558308
4473          1 2558740  5724264
4474          1 5724940  5749083
4475          1 5749226 12529544
>

https://www.dropbox.com/s/wei3i787c9nrz8q/WTSI-OESO_121_1pre.copynumber.caveman.txt.csv?dl=0
I have another file like below
              Chromosome         Position
rs62635286       1                 51479
rs75454623       1                 114930
rs806731         1                 30923

https://www.dropbox.com/s/tiqjhkrpfr9d9gu/WTSI-OESO_121_1pre.copynumber.txt?dl=0
In first  column I have some SNPs, how I can count the number of SNP in each range? For example based on the second table how many SNP are in range of 51479 to  817980?

Repost of https://www.biostars.org/p/376315/

We need to "merge with overlap/range" then "group by count", searching these terms should help you solve the problem.

Damian and Wouter have given you perfect solutions on BioStars. There is nothing more to add.

@FereshTeh I think people are getting a bit upset if you crosspost in different sites and it's not like you are a new user who just does not know. You are here for a while and you know how this site works. You must have noticed that crossposting is not very popular and if you do, at least you can explicitly mention it in the question.

Even better, if you get answers in one of the two places you ask, at least take the time to post them in the other. So if you get an answer here, also post it on your question at biostars and if you get an answer there, also post it here. That way, both communities benefit.

@FereshTeh I found a good article for you. Try always explain the problem (you do that moreless already), what have you tried to resolve it (this one you need to work on) and how to reproduce the problem you find (this one you do, but it would be nicer if you shared samples, not dropbox links).

There is like a thousand different ways how to achieve this. You could use a specialised software for this (like bedtools) or calculate it simply in R.
R solution: You can make a function that calculates the number of SNPs in a range and than you can apply it on all the ranges in the table with genomic ranges.
snp_table <- read.table('WTSI-OESO_121_1pre.copynumber.txt')
genomic_ranges <- read.table('WTSI-OESO_121_1pre.copynumber.caveman.txt.csv', sep = ',', col.names = c('Chromosome', 'Start', 'End'), stringsAsFactors = F)

how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(snp_table$Chromosome == coords[1] & snp_table$Position > as.numeric(coords[2]) & snp_table$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_snps <- apply(genomic_ranges, 1, how_many_in_range)

Now the table genomic_ranges has a variable that sais how many SNPs are in the range.

You know, the problem others have with your question is not about your coding skills, but more about the crossposting that really annoys the heck out of everybody, also your questions do not clearly show much research effort. It might be that you spent days trying to figure it out, but from your question, we can not tell, because you did not write what you have tried. Finally, I don't think self-hatred is useful to improve your coding skills, reading manuals and trying it out is way more efficient.


text:
Loading local FASTA file in igv.js
I'm trying to use igv.js for a very simple visualization of some genomic alignments, but I can't get past the first step.  I have a FASTA file that I want to use as the reference, but igv.js consistently gives an error when I try to load that file.  The file is a valid FASTA (have tested others from other sources to be sure - same error), and I can specify a remote https:// file, but not a local file.
Here's the HTML page I'm working on:
<!DOCTYPE html>

    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no">

    <title>Local FASTA Test</title>

Local File Test

and here's the FASTA file
>HSBGPG Human gene for bone gla protein (BGP)
GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT
ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC
ATGAGAGCCCTCACACTCCTCGCCCTATTGGCCCTGGCCGCACTTTGCATCGCTGGCCAGGCAGGTGAGTGCCCC
CACCTCCCCTCAGGCCGCATTGCAGTGGGGGCTGAGAGGAGGAAGCACCATGGCCCACCTCTTCTCACCCCTTTG
GCTGGCAGTCCCTTTGCAGTCTAACCACCTTGTTGCAGGCTCAATCCATTTGCCCCAGCTCTGCCCTTGCAGAGG
GAGAGGAGGGAAGAGCAAGCTGCCCGAGACGCAGGGGAAGGAGGATGAGGGCCCTGGGGATGAGCTGGGGTGAAC
CAGGCTCCCTTTCCTTTGCAGGTGCGAAGCCCAGCGGTGCAGAGTCCAGCAAAGGTGCAGGTATGAGGATGGACC
TGATGGGTTCCTGGACCCTCCCCTCTCACCCTGGTCCCTCAGTCTCATTCCCCCACTCCTGCCACCTCCTGTCTG
GCCATCAGGAAGGCCAGCCTGCTCCCCACCTGATCCTCCCAAACCCAGAGCCACCTGATGCCTGCCCCTCTGCTC
CACAGCCTTTGTGTCCAAGCAGGAGGGCAGCGAGGTAGTGAAGAGACCCAGGCGCTACCTGTATCAATGGCTGGG
GTGAGAGAAAAGGCAGAGCTGGGCCAAGGCCCTGCCTCTCCGGGATGGTCTGTGGGGGAGCTGCAGCAGGGAGTG
GCCTCTCTGGGTTGTGGTGGGGGTACAGGCAGCCTGCCCTGGTGGGCACCCTGGAGCCCCATGTGTAGGGAGAGG
AGGGATGGGCATTTTGCACGGGGGCTGATGCCACCACGTCGGGTGTCTCAGAGCCCCAGTCCCCTACCCGGATCC
CCTGGAGCCCAGGAGGGAGGTGTGTGAGCTCAATCCGGACTGTGACGAGTTGGCTGACCACATCGGCTTTCAGGA
GGCCTATCGGCGCTTCTACGGCCCGGTCTAGGGTGTCGCTCTGCTGGCCTGGCCGGCAACCCCAGTTCTGCTCCT
CTCCAGGCACCCTTCTTTCCTCTTCCCCTTGCCCTTGCCCTGACCTCCCAGCCCTATGGATGTGGGGTCCCCATC
ATCCCAGCTGCTCCCAAATAAACTCCAGAAG

Welcome to the site.  Have you reached to the IGV.js developers?  What is the "error" that you get

I have not reached out to the developers - I want sure where to do that.  The error says "cannot find chromosomename of undefined" or something similar.  The key being "undefined"which I interpret as, it did not load the specified file.

The exact error I get is a dialog box stating:  Cannot read property 'getHomeChromosomeName' of undefined

It looks as though the script is looking for a URL (fastaURL) so you could try specifying the file as file://directory/test.fasta
By default though, javascript cannot read local files from disk for security reasons. Are you sure it is even possible to specify a local file in your script in this way?  If yes, then the problem is likely that you need to specify the path correctly.

Thanks for the comment.  I tried specifying the file as you suggested, but no luck.  Perhaps it is the JavaScript limitation.  I suppose I respect why that design decision is made, but for in internal use only purpose, it seems limiting.

One way to get round it is to host the web page and the FASTA file together in your own local webserver directory.  OSX has a built in local webserver that you can switch on without too much difficulty - then open the page from your browser using http://localhost/test.html   You would then be able to refer to the FASTA file as test.fasta as your script does now (URL relative to the calling page).

That did the trick.  Both files must be on a server and it does indeed work.

just stumbled on this.  I'm the developer of igv.js.  Please feel free to open issues and questions on our git forum,  I don't often have time to browse other forums.
Browsers will not allow loading of local files by file path in javascript, if they did it would be a gigantic security hole.  This restriction includes "file://" urls.  Local files can be loaded by user action (clicking a button),  that is supported in igv.js but there is no button for it as igv.js is designed as a component to embed in pages, not an application by itself.  For an example see igv.org/app (https://github.com/igvteam/igv-webapp),  a minimal application wrapping igv.js.  You'll find "local file" options for genomes and tracks there.


text:
No refseq transcripts in 2nd half of Y chromosome for HG38
Obviously there is a gap in my understanding of how sex chromosomes are annotated. I've been working on some CNV calls and noticed that most of the 2nd half of ChrY is missing annotations and no reads are being mapped there. Yq12 seems to be all NNNNNNN in the hg38 reference.
Anyone know why?

Nevermind, Answered my own question. It's a repetitive region that isn't very well understood

The human Y chromosome, because it is enriched in repetitive DNA, has been very intractable to genetic and molecular analyses. There is no previous evidence for developmental stage- and testis-specific transcription from the male-specific region of the Y (MSY). Here, we present evidence for the first time for a developmental stage- and testis-specific transcription from MSY distal heterochromatic block. We isolated two novel RNAs, which localize to Yq12 in multiple copies, show testis-specific expression, and lack active X-homologs. Experimental evidence shows that one of the above Yq12 noncoding RNAs (ncRNAs) trans-splices with CDC2L2 mRNA from chromosome 1p36.3 locus to generate a testis-specific chimeric β sv13 isoform. This 67-nt 5′UTR provided by the Yq12 transcript contains within it a Y box protein-binding CCAAT motif, indicating translational regulation of the β sv13 isoform in testis. This is also the first report of trans-splicing between a Y chromosomal and an autosomal transcript.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1832090/


text:
cooler out of memory
I ran into the following error when using cooler to load a contact pairs file:
ubuntu@ip-172-31-8-42:/data/hg-tmp$ FILE=GSE109344_monkey_fibro_allValidPairs.txt.gz; sudo docker run -v $(pwd):/data $TAG cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} /data/$FILE /data/${FILE}.cool;
Traceback (most recent call last):
  File "/usr/local/bin/cooler", line 11, in <module>
    sys.exit(cli())
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 764, in __call__
    return self.main(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 717, in main
    rv = self.invoke(ctx)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 956, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 555, in invoke
    return callback(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/cooler/cli/cload.py", line 477, in pairs
    map(pipeline, reader),
  File "cytoolz/functoolz.pyx", line 493, in cytoolz.functoolz.Compose.__call__
  File "/usr/local/lib/python2.7/site-packages/cooler/create/_ingest.py", line 425, in _aggregate_records
    .aggregate(agg)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 1315, in aggregate
    return super(DataFrameGroupBy, self).aggregate(arg, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 186, in aggregate
    result, how = self._aggregate(arg, _level=_level, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 498, in _aggregate
    result = _agg(arg, _agg_1dim)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 449, in _agg
    result[fname] = func(fname, agg_how)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 432, in _agg_1dim
    return colg.aggregate(how, _level=(_level or 0) + 1)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 760, in aggregate
    return getattr(self, func_or_funcs)(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/groupby.py", line 1227, in size
    result = self.grouper.size()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 225, in size
    ids, _, ngroup = self.group_info
  File "pandas/_libs/properties.pyx", line 36, in pandas._libs.properties.CachedProperty.__get__
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 252, in group_info
    comp_ids, obs_group_ids = self._get_compressed_labels()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 272, in _get_compressed_labels
    return compress_group_index(group_index, sort=self.sort)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/sorting.py", line 379, in compress_group_index
    comp_ids, obs_group_ids = table.get_labels_groupby(group_index)
  File "pandas/_libs/hashtable_class_helper.pxi", line 1255, in pandas._libs.hashtable.Int64HashTable.get_labels_groupby
  File "pandas/_libs/hashtable_class_helper.pxi", line 199, in pandas._libs.hashtable.Int64Vector.resize
MemoryError: cannot allocate memory for array

The machine has 8Gb of RAM. Any ideas?

It seems to be some kind of a memory leak in pandas's read_csv, maybe due to it doing gzip-decompression.
Doing the decompression in a separate process and piping the result into cooler cload pairs should do the trick.
gzip -dc GSE109344_monkey_fibro_allValidPairs.txt.gz | cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} - /data/${FILE}.cool


text:
UCSC packages/tools ucsc-netchainsubset and ucsc-fatotwobit not working on GNU/Linux
I have installed some ucsc tools on my Mac computer and it worked fine:
Installation
conda create --name ucsc-netchainsubset377 --channel bioconda --yes ucsc-netchainsubset=377

Conda environment activation
source activate ucsc-netchainsubset377

Checking the tool's help
    netchainsubset
netChainSubset - Create chain file with subset of chains that appear in the net
usage:
   netChainSubset in.net in.chain out.chain
options:
   -gapOut=gap.tab - Output gap sizes to file
   -type=XXX - Restrict output to particular type in net file
   -splitOnInsert - Split chain when get an insertion of another chain
   -wholeChains - Write entire chain references by net, don't split
    when a high-level net is encoundered.  This is useful when nets
    have been filtered.
   -skipMissing - skip chains that are not found instead of generating
    an error.  Useful if chains have been filtered.

The same commands on GNU/Linux produced the following error:
netchainsubset
netchainsubset: command not found

even though the tool seems to be installed successfully:
Downloading and Extracting Packages
ucsc-netchainsubset- | 310 KB    | ############################################################################ | 100%
mysql-connector-c-6. | 4.4 MB    | ############################################################################ | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

The same happened with ucsc-fatotwobit and ucsc-liftup and I would not be surprised if this behaviour affected the whole suite of ucsc programmes installed through conda.
Anyone can suggest a solution?
Update 5 minutes after posting
Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.
ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.
Cheers terdon.
Cheers.

Please [edit] your question and tell us how you installed the tools. It sounds like they're simply not in your $PATH, and you need to either add the directory where you installed them to your $PATH or just call the program using its full path (/path/to/ucsc/netchainsubset).

Short answer: respect the case of the tool name
Longer story below
Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.
ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt
Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.

Yes, the default file system macOS uses, HFS+, is case insensitive. One of the very few file systems that are. That's why the case didn't matter on your mac. However, most file systems, including those used by Linux, are case sensitive and A is not the same as a.


text:
repeatscout freq table
I'm running Repeatscout. I built the l-mer table called myfile.freq of myfile.fa
Can anyone tell me what do they mean the second and third columns produced as output?
here I report an example:
AAAAAAAAGCGGGA  3       107776875
AAAAAAACTGTATG  10  83440519
AAAAAAAAGGCGTA  3   41037187
AAAAAAACTTGAAT  7   94493612
CATACATGCATGCA  1065    125671338
CATACATGCTTGAA  7   121799834
AAAAAAATCATGCA  10  95493021
AAAAAAAGTCCAGT  3   125127980
AATTCACATGTATG  7   102505668

I read the tutorial that explains how to run the program: https://bix.ucsd.edu/repeatscout/readme.1.0.3.txt.
It says: "First, build_lmer_table creates a file that tabulates the frequency of all l-mers in the sequence to be analyzed"

Welcome to the site Jonny. Did you read the documentation of the software you are using, or any tutorial where it is used? Did you try searching on the internet for this format/output?

Yes, I searched and I read some manuals but I really don't understand what they mean. I know that it's built this table of l-mer from the original genome fasta file

If you could link and cite what parts you don't understand we could help you better

" First, build_lmer_table
creates a file that tabulates the frequency of all l-mers in the
sequence to be analyzed."  I SUPPOSED that the second column could indicate the number of l-mer, but I don't understand what the third column is showing

Where are you taking this from? Could you provide also the link?  Do not forget to [edit] the question to include this relevant information


text:
Calculating p-value for introns which are retained (or not) between 2 conditions
I am trying to calculate p-values and ultimately FDR values for RI introns for 2 conditions of experiment I have(control and knockdown)
For every condition I have 3 biological replicates each: 3 control replicates and 3 knockdown replicates.
This is the header line of my input file.
chrom   Istart  Iend    strand  avg-L_R_counts-1    Sk_counts-1 avg-L_R_counts-2    Sk_counts-2

-1 and -2 are 2 sample conditions control and knockout respectively
avg-L_R_counts = for a particular intron, this is the average of reads 5bp up stream and downstream of the 5' end and at the 3'end of the intron
Sk_counts = This is the number of reads spanning exon-exon junction for that intron.
I have approximately 350 introns(this is yeast). I can get the Percent Spliced In values but I am interested in getting a p-value and then I can calculate the FDR values.
Can anyone point me as to what approach I should take to calculate a p-value, or what test should I use?

What is your null hypothesis? What do you want to test? That one condition has more introns than the other?

My null hypothesis is that there is no intron retention in either conditions

Here is how we do this:
First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.
transcript 1 |>>>>>>>>>|----------|>>>>>>>>|---|>>>>|
transcript 2 |>>>>>>|-------------|>>>>>>>>|---|>>>>|
transcript 3 |>>>>>>|-------------|>>>>>>>>>>>>>>>>>|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|

this can be achieved using cgat gtf2gtf --method=genes-to-unique-chunks from the cgat-apps toolkit if you don't want to code it yourself.
I call
chunks 1, 4 and 6 constitutive exon
chunk  3 constitutive intron
chunk  2 alternate intron
chunk  5 annotated retained intron.
Now use featureCounts to count the number of reads in each chunk.
Test the chunk counts for usage using DEXSeq.
Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.

If you have your transcripts in BED format, you can use bedops --partition to make chunks.

Did you look into literature? I am not an expert in this field, but a small search on google leads me to IRFinder (ref.), a tool that can assess intron retention. My advice is to see if this tool can help you, before reinventing the wheel...

You might take a look at the MISO algorithm for inspiration, or as a method for calculating confidence intervals of alternatively spliced isoforms: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3037023/


text:
What is the `table` parameter in the UCSU ldHgGene tool?
I want to lift over the annotation from E. coli release 29 to E. coli release 42. I have followed the steps for an annotation liftover mentioned here. I got the chain file for the final step and ran:
liftOver -gff Escherichia_coli.GCA_000967155.1.gtf Chromosome.chain conversions.gtf unmapped.gtf
However, I get the following warning:
WARNING: -gff is not recommended.
Use 'ldHgGene -out=<file.gp>' and then 'liftOver -genePred <file.gp>'

I would like to heed the warning. However, I don't know how to use the ldHgGene tool. Here is its help page:
ldHgGene - load database with gene predictions from a gff file.
usage:
     ldHgGene database table file(s).gff
options:
     -bin         Add bin column (now the default)
     -nobin       don't add binning (you probably don't want this)
     -exon=type   Sets type field for exons to specific value
     -oldTable    Don't overwrite what's already in table
     -noncoding   Forces whole prediction to be UTR
     -gtf         input is GTF, stop codon is not in CDS
     -predTab     input is already in genePredTab format
     -requireCDS  discard genes that don't have CDS annotation
     -out=gpfile  write output, in genePred format, instead of loading
                  table. Database is ignored.
     -genePredExt create a extended genePred, including frame
                  information and gene name
     -impliedStopAfterCds - implied stop codon in GFF/GTF after CDS

From what I understand, I am trying to give the liftOver tool the annotation in GenePred format instead of GTF format. So I need to convert the GTF format into a GenePred format and that's what the ldHgGene tool is for. So the command should be something like:
ldHgGene Escherichia_coli.GCA_000967155.1.gtf -out=Escherichia_coli.GCA_000967155.1.gp

but this returns nothing. The -out option does mean that the database parameter is not needed but the tool still seems to need the parameter table.
What is the table parameter?

I stumbled across this issue after getting the same warning message from liftOver. By using the -out option of lgHgGene, the database parameter is ignored but must still be provided. The table parameter should be the second column from your GTF file. Note that you must also specify the -gtf option (otherwise it expects GFF format). In my case, the GTF file looks like:
chr1    mm10_knownGene  exon    3073253 3074322 0.000000        +       .       gene_id "ENSMUST00000193812.1"; transcript_id "ENSMUST00000193812.1";
chr1    mm10_knownGene  exon    3102016 3102125 0.000000        +       .       gene_id "ENSMUST00000082908.1"; transcript_id "ENSMUST00000082908.1";
chr1    mm10_knownGene  exon    3206523 3207317 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3213439 3215632 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3205901 3207317 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
chr1    mm10_knownGene  exon    3213609 3216344 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
...

By default, ldHgGene outputs a leading 'bin' column which won't work with liftOver; this can be avoided with the -nobin option. So my command was:
ldHgGene -gtf -nobin -out=mm10.GENCODE_VM20.gp ignored mm10_knownGene mm10.GENCODE_VM20.gtf

The resulting genePred output can be fed directly into liftOver -genePred:
ENSMUST00000193812.1    chr1    +       3073252 3074322 3074322 3074322 1       3073252,        3074322,
ENSMUST00000082908.1    chr1    +       3102015 3102125 3102125 3102125 1       3102015,        3102125,
ENSMUST00000162897.1    chr1    -       3205900 3216344 3216344 3216344 2       3205900,3213608,        3207317,3216344,
ENSMUST00000159265.1    chr1    -       3206522 3215632 3215632 3215632 2       3206522,3213438,        3207317,3215632,


text:
CDS length for each human gene
Does anyone know where and how could I download a list of all human genes and the length of the coding sequence for each gene? Is it possible to do this on the NCBI site, ensembl?

Which coding sequence? I mean, do you just want whichever has been designated the 'canonical' transcript or do you want all possible isoforms?

Hi terdon, thanks for the quick reply! Yes exactly, the canonical transcript is good enough!

Michael G. suggests to take a look at relevant front-end, NCBI's eFetch. Which is supposedly perfect for what you need.

While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them).
Essentially, you just need to go to BioMart, and

select "Ensembl Genes 96" (the number will change if the version changes) as the database and "uman Genes" as the dataset.
Click on "Filters", and set Gene type to coding and Transcript type to protein_coding.
From "Attributes", select whatever you want to see. The "CDS Length" is under "Structures".

Ensembl has an FTP site that allows you to select and download only the coding sequences from many different genomes.
https://useast.ensembl.org/info/data/ftp/index.html
To determine the length of those sequences, download the associated gtf or gff3 annotation file. The annotation file is tab delim. The fourth and fifth column represent genomic loci of the annotated region. Subtract the amount in the fourth column from the amount in the fifth column to yield the length of all the annotated features.
You can easily do this in R after loading the file in the environment using the Magrittr library. The following code will create a new column called gene_length with the associated gene lengths.
install.packages("magrittr")
# this only needs to be done once
library(magrittr)
# must be run each time the library is neaded
annotation.gtf <- read.table("path/to/annotation.gtf")
annotation.gtf$start <- annotation.gtf[,4]
annotation.gtf$end <- annotation.gtf[,5]
annotation.new-column.gtf <- annotation.gtf %>%
          mutate(gene_length=end-start)

Go to the NCBI homepage and search for 'human genome' (the dropdown menu on the left-hand side of the search box should be 'All Databases'). In the result page, click on the 'Download' button; choose 'RefSeq' as the 'Source database' and 'Feature table' as the 'File type' as shown in the screenshot below:

This will download a genome_assemblies.tar file which contains the feature table file, GCF_000001405.39_GRCh38.p13_feature_table.txt.gz. This file has all of the annotated features in a tab-delimited format with the following fields:
              # feature [  1]: gene
                  class [  2]: transcribed_pseudogene
               assembly [  3]: GCF_000001405.39
          assembly_unit [  4]: Primary Assembly
               seq_type [  5]: chromosome
             chromosome [  6]: 1
      genomic_accession [  7]: NC_000001.11
                  start [  8]: 11874
                    end [  9]: 14409
                 strand [ 10]: +
      product_accession [ 11]:
   non-redundant_refseq [ 12]:
      related_accession [ 13]:
                   name [ 14]: DEAD/H-box helicase 11 like 1
                 symbol [ 15]: DDX11L1
                 GeneID [ 16]: 100287102
              locus_tag [ 17]:
feature_interval_length [ 18]: 2536
         product_length [ 19]:
             attributes [ 20]: pseudo

You can parse this file to obtain the CDS lengths of all annotated proteins. As an example, you can get a list of all annotated RefSeq proteins and their lengths as follows:
zcat GCF_000001405.39_GRCh38.p13_feature_table.txt.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($1~/^#/ || $1=="CDS")' \
    | cut -f6-11,13,15,16,18,19

The output of the above command will have the following fields:
             chromosome [  1]: 1
      genomic_accession [  2]: NC_000001.11
                  start [  3]: 69091
                    end [  4]: 70008
                 strand [  5]: +
      product_accession [  6]: NP_001005484.1
      related_accession [  7]: NM_001005484.1
                 symbol [  8]: OR4F5
                 GeneID [  9]: 79501
feature_interval_length [ 10]: 918
         product_length [ 11]: 305


text:
How to get Uniprot and swissprot cross references using Biopython
I am struggling to use biopython to gather cross references such as GO annotation from Uniprot/swissprot text files ("DR" in the .txt files) using biopython.
From the source code:
def _read_dr(record, value):
    cols = value.rstrip(".").split('; ')
    record.cross_references.append(tuple(cols))

This leads me to think that the following code should work, as this is how other features like annotation and features work.
for record in SeqIO.parse(filename, "swiss"):
    print(record.cross_references) #does not work, should return a list of DR rows
    print(record.features) #works, returns a list of parsed FT rows

I believe the attribute you should be looking for is dbxrefs:
[x for x in record.dbxrefs if x.startswith('GO:')]

Apparently SeqIO handles uniprot text files differently than the SwissProt module:
from Bio import SwissProt
record = SwissProt.read(open('single_record_from_uniprot.txt'))
record.cross_references

So you should change you code to use the SwissProt module to get the behavior you expect


text:
Repeatscout installation problem
I've installed Repeatscout from https://bix.ucsd.edu/repeatscout/.
I've also installed TRF and nseg required by the program.
I followed the steps:
build_lmer_table -l 16 -sequence ~/file_fasta/od_1000.fa -freq od_1000.freq

 RepeatScout -sequence od_1000.fa -output od_1000_repeats.fa -freq od_1000.freq -l 16

and the I put this command:
cat od_1000_repeats.fa| filter-stage-1.prl >od_1000_repeats_filt1.fa

but this message appears:
bash: filter-stage-1.prl: command not found...

Why? It is all correctly installed, I can't understand the problem

If you're in the RepeatScout directory you need to execute the Perl script with a ./ prefix, like so.
cat od_1000_repeats.fa | ./filter-stage-1.prl > od_1000_repeats_filt1.fa

If you're not in the RepeatScout directory, you need to provide the full relative or absolute path when executing the script. Alternatively, you can put the RepeatScout directory in your $PATH environmental variable.
These are all standard UNIX/Linux concepts that many bioinformaticians take for granted, for better and for worse. If you're lost in all the terminology, a crash course in UNIX basics is probably where you would want to start.

I installed repeatscout on my PC, I created the Repeatscout directory and It worked without errors. The problem is that it can't execute a file larger more than 800 mb in a 16 gb ram of my pc. So I had to use the cluster but I don't know where is the Repeatscout directory.(If I put "locate Repeatscout" I don't find it). I would to find the path to check if all the programs are correctly installed..

It sounds like you need to get some help from your system administrator.


text:
What does PCA mean on GWAS
I understand what GWAS is and I'm able to perform certain tests with the p-values, etc. But what I am having a hard time wrapping my head around is what PCA on GWAS means.
So let's say I have 100,000 individuals and genotype data for 10 million SNPs. And I'm analyzing this for 10 different diseases. I also have p-value information for each of the 10 diseases for the 10 million SNPs.
At my work, they were throwing around the following phrases. I'm not quite sure what exactly it means.

"Perform PCA on GWAS." Does this mean, I'm performing PCA on the genotype data or p-value data? I've done PCA on classic machine learning problems like "iris." But I'm not able to quite wrap my head around what PCA information tells me in terms of this genetic data.
"PCA on individual diseases vs. PCA on entire GWAS" What exactly is the difference. I realize this is a broken statement, but they said we have two sets of PCs. One on the individual diseases which are unique, and then PCs on the entire GWAS for all 10 diseases.

Sorry if my questions is not properly asked, I realize it's probably heavily abusing notation.
Added clarifications:

So, my matrix has 100k rows and 10 million columns. Then, PC1, for example, is a vector with 100k values, each value corresponding to each individual. What exactly does PC1 value that's associated with the first individual tell me with regards to the PC1 value that's associated with 2nd individual, and so on. My current understanding if you think of a number line that represents PC1, then each of the PC values (corresponding to each individual) represents where on the number line that individual falls.
How exactly does this eliminate/account for variation in population stratification?

The questions are quite easy. What's your sample distribution? What is the sample size of healthy individuals and diseased individuals per disease?

I have about 95,000 healthy individuals and 5k diseased. @MichaelG.

Oh.. hugely imbalanced, that is difficult. You mention 10 diseases, is that 500 per disease? This last question is quite important

Sorry, I should have been more clear. For each of the 10 diseases, we have event data for 100k individuals. So for the first disease, 95k of them are healthy. But for instance, one of the traits has 30k individuals with the trait and 70k without. And we're looking at the same 100k individuals across all 10 diseases. @MichaelG.

I think I get it. Let compose an answer. The overlapping populations sound a bit tricky

From my memory of what a statistician told me, a PCA aims to determine independent linear combinations of variables (i.e. genotypes) that account for the most variation in the dataset. With 10 million SNPs, the vectors describing the genotypes of all individuals reside somewhere inside a 10-million-dimension hypersphere, and the first principle component determines the best vantage point / angle from outside that sphere that will create the greatest spread of data when data is projected onto a perpendicular / orthogonal axis. The second principle component is the perpendicular / orthogonal axis that will create the second-greatest spread of data.
Aside: There's some wonderful linear algebra (which I've forgotten) that makes it very quick (computationally) to do this angle calculation, exploiting the nature of linear transformations and orthogonality.
The point of the PCA for GWAS is to identify (and preferably exclude) differences that might be due to the structure of populations, rather than the disease of interest. It's fairly common for researchers to do a PCA on the entire set of genotype data, then use the first 5, 10, or 20 principal components as covariates in the association model.
There's an underlying assumption that the majority of genetic variation is due to population structure - if this is not the case, then the disease would have to be extremely omnigenetic, with almost the entirety of the genome contributing to the disease. The approach also won't deal with all the population structure (because there may be non-linear contributions of genotypes to population structure), but it will get rid of the largest problems.
Additional to this, I'd recommend the following things:
Structure Visualisation
Carry out a STRUCTURE analysis or PCA on the combined cases and controls, and display the results of the first two principle components, labelling, colouring, or otherwise indicating both groups. As an additional check, it may be useful to do addititional PCAs for cases and controls separately.
If they appear different, it suggests that the cases and controls were drawn from separate population clusters, which should be commented on in the discussion of the paper. There's no need to panic and find more data, just make the observation.
Subsampling
Carry out association tests on an equal-sized subsample of both datasets. This will mean that area-under-the-curve statistics can be properly interpreted. If you've got unbalanced datasets, use a little bit less than the smaller of the two groups (e.g. for 5k diseased and 95k healthy, use a random subsample of 4.5k from both groups).
Multiple testing / Internal validation
Run the association tests multiple times (with different subsamples), recording the individuals that were subsampled for future reference / reproducibility.
If an association is only present in one of the many different subsamples, it should be ignored as a spurious result. It means that there is some association that is not generalisable to all (or most) of the cases and controls. These variants might be interesting in themselves on an individual level, but not in a general population.
I explain some of this in more detail in the following preprint (which I'm not planning to submit to any journal, due to it being a bit stale and needing a lot of work to fix up):
Bootstrap Distillation: Non-parametric Internal Validation of GWAS Results by Subgroup Resampling

I've added in a clarification, please see OP

Note that in this context "independent" should be "orthogonal". There's another dimensionality reduction method called Independent Component Analysis where the variables are linear combinations of input variables that are maximally statistically independent of each other.
Everything else here is good advice

I did mention orthogonal, just not in the introductory sentence.

PCA = principle component analysis and a multivariate statistic, today it is trendily retermed "unsupervised learning" and here is likely being deployed for individuals within your data set. It works by identifying the maximum variance within multidimensional space, shearing it and describing this as the first principle component. The second principle component is sections the variance at 90o to the first in all axes (orthogonal is the posh word), and so on ... The maximum information is within the first two components and they are sometimes plotted on x and y axes to look for potential clusters within your data. In this case the first few components are all that is needed because you not identifying the genes of interest just checking the data looks sound.
Point 1
PCA on all GWAS is attempting to assess that the genetic variation between diseased and healthy people is evening spread across the genome. For example if your diseased population came from e.g. the indigenous people of Western Samoa (purely as an example), which has a very defined genetic heritage (at least for HLA) then most of the SNPs will be unrelated to the disease your are attempting to identify. This will mean you have to redefine the healthy population, the null distribution, to screen out the artefacts. In this case a better representation would be to select healthy individuals from Western Samoa as your comparison. So within the first two components you want to see lots of mixing between your patients vs healthy people. If the patients are clumping (clustering) thats bad here and you have to find out why.
PCA has weaknesses BTW for complicated reasons which are boring, but it is often not assessed directly but via a tSNE analysis and then plotted on two axes.
Point 2
If PCA from point 1 looks cool, then you might split the data set into diseased and another one of healthy and perform a PCA on each but now look at the SNPs ... not the people. I hope that makes sense. If you compare the PCA plot of SNPs between these two groups, the differences MIGHT and I stress might highlight the SNPs of interest, i.e. the SNPs causing your disease. However and its a big however the number of SNPs you are dealing with is huge and if you are manually going through all that data to assess anomalous clusters, or outliers to make a list of candidate SNPs for further analysis (@gringer describes some of the further options above) -  I wish you good health and long life, and you will need it because you will be very old at the end of that optical analysis.
What you really do is an additional analysis (I would deploy tSNE) to assess the differences between the two PCA plots and end up with a third plot indicating the differences between all the other PCA plots. In this context it is complicated to deploy and would be nested within a tSNE analyse of all components which is what I alluded to in point 1.
Point 2 - issues Generally if someone thought I was going to eye-ball 50000 SNPs between two graphs, or throw them back a reanalysis which conceptually made their brain hurt ... their gonna have conceptualise real hard.
Summary Having said all that PCA is merely an exploratory tool and is not the defining significance test required, to find the "one SNP" which will explain your disease.

The imbalance of the data set makes advanced SNP identification statistically difficult, but the 30:70 split is okay. However a subsample of the healthy population, which is what Point 1 is all about, would likely correct for this imbalance.

On the question of mixing Point 1, what you don't want to see is your diseased patients forming a "blob", or even a circle, or any sort of defined cluster against the spread of healthy patients on a scatterplot of PC1 Vs PC2. Perhaps I should have said if your patients don't cluster on the output to Point 1 that's good. There are statistics that identify clusters on a PC1 PC2 scatter plot, and you want it to fail that test.

I'm a bit confused on point 1. I've added clarification in my question.

What exactly would "a lot of mixing between patients and healthy people" look like in terms of PC plots? How would I identify this?

Nevermind, I totally understand what you're saying now, thank you!

I'll add a more detailed explanation and think of an example


text:
PAML branch-site model: Identical lnL value between null and alternative models?
I have been running PAML, via codeml, to try to identify positive selection in certain branches of my tree for one gene (but will eventually be trying to automate the analysis over multiple genes). This was done by setting model = 2 and NSsites = 2 in the codeml.ctl file. I attempted to run both the null and alternative models by setting fix_omega = 1 in the null.ctl and fix_omega = 0 in the alternative.ctl.
However, after running both models, they have returned exactly identical log-likelihood (lnL) values. I wasn't expecting to find anything significant, but the fact that they were identical struck me as odd. Out of curiosity, I altered the branch of interest in the newick tree file and reattempted the analysis, but again, the lnL values remain the same between the null and alternative models (although slightly different to the value returned with the different branch selected).
I have since repeated this many times, choosing a different branch(es) each time to see if I can get a difference between the null and alternative, but no luck. Am I missing something obvious? Should I really be getting the same exact result between null and alternative models, regardless of where in the tree I place #1? To reiterate, I wasn't expecting to find significant results, but was at least expecting (non-significantly) different numbers for lnL...

The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.
Possible explanation fix_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.
Goto debugging As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,
time codeml ./file

The thing to remember is fix_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.
Nitty gritty debugging
1 Please check you have
seqtype = 3

If not, then that is your answer, omega was never being set.
2 For debugging only make sure you have,
model = 1

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below
3, Then try a different NSsites model,
NSsites = 8

This should produce real differences. The NSsites options are,

0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w9:betaγ
10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1; * 13:3normal>0

w here is for omega
4 The proof of your system is that when you
fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

This should produce likelihood variation

Finally
model = 2

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it.  I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.

Thank you Michael, you've given me something to try. If I'm still struggling after debugging, I'll update this question.

Any updates @PollardMD ?

Apologies for the late reply. I managed to get it working eventually. It looks like some genes/branches that I am looking at truly do produce identical values for both null and alternative models. When looking at other genes of interest, I have been able to get non-identical results from the null and alternative models. I guess I just wasn't expecting the lnL values to be perfectly identical, and that threw me off.


text:
GSEA - Gene Set Enrichment Analysis - in R studio
I have been trying to conduct a classic GSEA using the broad institute software and can not figure out how to compile it in R for the life of me. I am aware that there exists other packages in R that are almost the same, such as fgsea however, I would like to use this one first and then compare it with these other methods.
I have already generated a ranked list for my contrast of interest and save it as a .rnk file. Here is what I tried to do:
gseaCallH <-
  paste("java -Xmx2G -cp C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7",sep="")

testGSEA <- system(command = gseaCallH, intern = TRUE, ignore.stderr = TRUE)

When I run the system command I get the following error:
had status 1
At first I thought maybe it had something to do with JAVA not being properly installed however, I have installed it and still nothing. Maybe the installation was not done correctly? I have tried running it directly from cmd (I'm running windows 10) and got the following error:
Error occurred during initialization of VM
Could not reserve enough space for 2097152KB object heap

Can anyone advise on how to proceed?

Are you sure you have the specs to run the application? Did you check the Java version?

BTW, do you have the rJava package installed? It requires to configure R to be able to use Java (basically running R CMD javareconf) maybe it helps too. Because the way you run it now from R is sending it back to the command line, so you should get the error from the command line

If by specs you mean the software then yes, I have downloaded the desktop version of the application as well

Yes, then it might be a problem for their support system... Good luck

I have answered my own question however, it is specific to windows users.. If you get the same error from the cmd as I did, you must manually put in environment variables related to Java. Do as the user Mohamed did in this SO post: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for
Following this, you must avoid using paths such as C:/absolute_path and instead use /absolute_path, if you do the former, you will get the exception: java.net.MalformedURLException: unknown protocol: c
In the end, after manually entering the environment variables and running the following command: java -Xmx2G -cp /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file /Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx /Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out /Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7 it has worked.
Where the command can be executed directly from the command line or in R using the system function.
Note: A '#' separated by a space must be present at the beginning of your .rnk file for the program to work. Therefore, .rnk files must have the following format:
# Genes  Coeffecients
TSPAN6    0.2345676
etc. for every gene.

Are you calling java with sufficient memory? Try increasing -Xmx2G to -Xmx8G, provided your PC has 8 GB RAM.

Good call. I tried this and got the exception Invalid maximum heap size: -Xmx8G

a quick google search suggests this might be due to 32 bit java vs 64 bit java?

So it actually had to do with the environment variables for Java in my computer. I had to manually insert some of them - _JAVA_OPTIONS and JAVA_HOME. Here is a SO post that was helpful: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for
Unfortunately I'm still getting an error - a different one though - malformed url exception, unknown protocol C. So I'm getting closer but still no luck..


text:
Generate table for total number of SV events per sample
I have called and annotated structural variants and I have a table by name of samples and structural events like below
Sample  Type
t_005   DEL
t_005   INS
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   BND
t_005   DEL
t_005   DEL
t_005   BND
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL

How can I count the number of events to get the below table?
> head(SV_Counts[,1:4])
               t_005      o_005       t_013         o_013
    DEL          27          40          25          30
    DUP          11          10          13          22
    BND          23          20          10          12
    TRA           3           2           4           4
    >

I tried this from google but I failed to adapt my data to that
    > melt(table(a))
   Sample Type value
1   o_036  BND   280
2   o_005  BND     8
3   o_013  BND   160
4   o_021  BND   144
5   o_036  BND     0
6   t_005  BND    15
7   t_013  BND   108
8   o_036  DEL    63
9   o_005  DEL    15
10  o_013  DEL   121
11  o_021  DEL   243
12  o_036  DEL     0
13  t_005  DEL    40
14  t_013  DEL    76
15  o_036  DUP    34
16  o_005  DUP     4
17  o_013  DUP    57
18  o_021  DUP    88
19  o_036  DUP     0
20  t_005  DUP    14
21  t_013  DUP    48
22  o_036  INS    58
23  o_005  INS     0
24  o_013  INS     4
25  o_021  INS    10
26  o_036  INS     1
27  t_005  INS     1
28  t_013  INS     2

But hereafter I don't know what to do

So you have data in "long format" and want to convert it to "wide format". You can use spread() from tidyr for that:
require(tidyr)
# Assuming your input is called `foo`
spread(foo, Sample, value)


text:
Frequency of specific viral sequence in .BAM or .fastq
I was wondering if anyone had any experience 'counting' the frequency of a particular viral sequence in an individual's fastq sequence. So basically, counting the number of occurrences of a particular sequence given someone's fastq file.
The thing is, we must speak with the sequencing company regarding what format we'd like the files and must choose between .bam and .fastq - both would be extra $. Usually I'd choose .bam out of convenience but I was wondering, given that it has already been aligned the human genome I might not be able to proceed with this 'counting' process. So, my question is: if I were to use samtools - or any other program - to convert the .bam files to .fastq, is there any possibility that I could lose such information - e.g. viral seqs - that would have otherwise been available should I have just gotten the .fastq files directly?

Rule of thumb is: never go with FASTQ if you can avoid it, it’s a net waste of space. Post-alignment BAMs are as small, or smaller, than gzipped FASTQ, despite containing more information.

Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.

Thank you !! Will check with the company to be sure

If you are interested in a particular viral sequence, it seems wasteful to request bams that are aligned to the human genome, but you could request bams aligned to the human genome with unaligned reads retained. As an alternative, you can generate unaligned bams, where all of the same reads are present that would be present in a fastq but in a bam format. You can revert this to a fastq or align the fastq directly to a viral reference yourself. Aligning to small references is fast.


text:
WGCNA error codes in Network Construction & moduel detection
I am using WGCNA and code has worked with no problems with some of my datasets however with most recent attempt when completing gene network & identification of modules. I get following error code...
net__30= blockwiseModules(t5, power=6,
                              TOMType = "signed", minModuleSize=30,
                              reassignThreshold =0, mergeCutHeight = 0.25,
                              numericLabels = TRUE, pamRespectsDendro = FALSE,
                              saveTOMs = TRUE,
                              saveTOMFileBase = "PBCreponseTOMB",
                              verbose= 3)

I come up with error code saying:
            Number of items to replace is not a multiple of replacement length.

Does anyone please know what this indicts to? I have looked on all the WGCNA tutorial pages but cannot find this error message anywhere.
Thank you kindly

Could be a bug in the code. You will need to provide more detail. When you get the error, type traceback() and post what you get. Also, post the output of sessionInfo().

Error code:     Error in clustDist[ii, ] <- if (intNetworkType == 1) 1 - abs(cr1) else 1 -  :
number of items to replace is not a multiple of replacement length

sessionInfo()    R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.5 LTS
Matrix products: default

3: sizeRestrictedClusterMerge(datExpr, clusters = membership, clusterSizes = clusterSizes, centers = centers, maxSize = preferredSize, networkType = networkType, verbose = verbose, indent = indent) 2: projectiveKMeans(datExpr, preferredSize = maxBlockSize, checkData = FALSE, sizePenaltyPower = blockSizePenaltyPower, nCenters = nPreclusteringCenters,  verbose = verbose - 2, indent = indent + 1)

1: blockwiseModules(B_CC_4, power = 6, TOMType = "signed", minModuleSize = 10,       reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "PBCreponseTOMB",
verbose = 3)

I can't add more of the sessionInfo() as too many characters.. can you let me know which parts are key to review. Thank you so much for your kind assistance. I am abit lost with the error and re-run over and over again and re-formatted my raw data, cant see any errors there, thank you

@holly You can [edit] your question to include these details and other that might be needed to give a good answer to your question (And help future readers to find it, as comments are not indexed by searcher)

Please try installing and using this development version of WGCNA: https://www.dropbox.com/s/fuu73li4ajjfr3z/WGCNA_1.67-90.tar.gz?dl=0 and see if it solves the issue. Thanks again for reporting the issue.

Thank you kindly for all your kind help, very much appreciate it

Looks like a bug in the code which I will try to hunt down and fix. In the meantime, I would suggest playing with maxBlockSize argument to blockwiseModules. Try increasing it as much as your available RAM allows (see the paragraph "A second word of caution concerning block size" on page 6 of WGCNA tutorial I, section 2c, for some guidance on how to set maxBlockSize). Hope this helps.

Thank you so much for looking over this for me, I very much appreciate your time and effort. I will try and play around with the maxBlockSize as suggested. Thank you kindly

I also got this error:
Did you manage to solve the issue @PeterLangfelder > > sessionInfo() R version 3.5.1 (2018-07-02) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 7 x64 (build > 7601) Service Pack 1 > > Matrix products: default > > locale: [1] LC_COLLATE=English_New Zealand.1252 LC_CTYPE=English_New > Zealand.1252 LC_MONETARY=English_New Zealand.1252 [4] LC_NUMERIC=C > LC_TIME=English_New Zealand.1252 > > attached base packages: [1] stats graphics grDevices utils > datasets methods base > > other attached packages: [1] WGCNA_1.67 fastcluster_1.1.25 > dynamicTreeCut_1.63-

Please try the newest version (1.68) of WGCNA from CRAN. This bug should be fixed.


text:
Finding matched genes for each genomic range
I have called structural variants and now I have
Sample  Type    Chrom1   Start1        End1 Strand1
t_005   DEL 1   28677074    28677079    +
t_005   INS 1   186383731   186383738   +
t_005   DEL 1   245902589   245902590   +
t_005   DEL 2   56345643    56345645    +
t_005   DEL 3   59766214    59766217    +
t_005   DEL 3   60270545    60270548    +

I have an annotation file like
> head(sv[,1:5])
  ENSG00000233440 13 23708313 23708703   HMGA1P6
1 ENSG00000207157 13 23726725 23726825    RNY3P4
2 ENSG00000229483 13 23743974 23744736 LINC00362
3 ENSG00000229483 13 23743974 23744736 LINC00362
4 ENSG00000252952 13 23791571 23791673  RNU6-58P
5 ENSG00000235205 13 23817659 23821323  TATDN2P3
6 ENSG00000232849 13 93708910 93710179 LINC00363
>

How I could relate each gene name from annotation file to my genomic range file?

EDITED

Sample  Type    Chrom1  Start1  End1    Strand1 Chrom2  Start2  End2    Strand2
t_005   DEL 1   28677074    28677079    +   1   28706324    28706329    -
t_005   INS 1   186383731   186383738   +   1   186383731   186383732   -
t_005   DEL 1   245902589   245902590   +   1   246007384   246007385   -
t_005   DEL 2   56345643    56345645    +   2   56345829    56345830    -
t_005   DEL 3   59766214    59766217    +   3   59785325    59785328    -
t_005   DEL 3   60270545    60270548    +   3   60705591    60705594    -
t_005   DUP 3   60744344    60744345    -   3   60806306    60806307    +
t_005   DUP 3   60882109    60882112    -   3   60934155    60934158    +
t_005   DEL 3   77973091    77973093    +   3   77991398    77991400    -
t_005   DEL 3   95972796    95972800    +   3   95974511    95974515    -

Hey Feresh, please use dput to present a reproducible example for both files then I can give it a shot. You just want to merge the sv file based on the range of start1 and end1 correct?

Yes likely I only need to know which gene locate in each genomic range. Genes are in sv file and my called structural variants are in the first file.

so why does end1 just have '+' and strand1 is empty?

Because structural variant data are too weird, I have added full called SV in edited post

I won't be able to help you without you posting a reproducible example using dput.. Please do this with a small snippet of your data as you have provided visually.

You received an answer from Kamil for a similar question here: https://bioinformatics.stackexchange.com/a/8532/650
You can adapt the code they gave you to compare Chrom1 == sv[2] & Start1 >= sv[3] & End1<=sv[4] to find events that happen within a gene. For events spanning genes, you're going to need to write a collapsing function for the gene names and a different way to compare overlaps instead of Start >= sv[3] & End <= sv[4].


text:
Merging transcriptomes coming from different experiments
I'm planning to build a transcriptome by pooling all existing transcriptomes in SRA for a non-model species (which has no reference genome) to study differentially expressed genes and the like. It might worth mentioning that the transcriptomes I want to pool mostly come from different tissues and sampling sites. These transcriptomes would come from five SRA BioProjects plus one project from my laboratoryfor a total of 102 single transcriptomes coming from about 20 individuals (plus three pooled larvae transcriptomes). I'm especially interested in differentially expressed genes in different tissues, on genes specific to different tissues and on the UTRs of differentially expressed genes.
My idea is to build a common reference transcriptome against which I'd re-assemble the individual transcriptomes, and only then I'd get differential expression data etc.
I'm looking for literature (or personal experience) dealing with assemblies of transcriptomes / comparing transcriptomes coming from different datasets in order to understand whether this is a good idea or not and get a grasp of what technical and biological problems I could face. I'm (vaguely) aware of the batch effect and of one way to remove it, but I'd like to read papers dealing with my same problem (different SRA databases, without reference genomes). Up until now I found just projects dealing with microarray data.

R's Bioconductor (obviously). There is a lot of experience on this board in exactly this area using NGS, albeit with the focus on "model" species (humans). No reference genome, sounds a bit scary.

I don't have personal experience, so I won't post an answer. Be aware the with ComBat you could introduce a new batch effect (it is not necessary to remove them in order to work with them, you can take it into account). The other problem is that you will have several conditions (even the same tissue is not the same at the first year of the organism that just a year before dying).

However the most important thing is what do you expect to learn? The differentially expressed genes don't mean anything if you don't know what you are comparing or if you are not looking for something very specific (healthy vs disease, one tissue vs the other...). Please [edit] the question to explain what is your biological question and it will be easier to help. (Also if you provide the number of projects of your organism and the number of total samples it would be helpful to us too)

@llrs I edited the question with the things you asked and I clarified it a little bit (I was too much into my project to notice that I had skipped some fundamental information...). I'm aware I will have different conditions, that's why I'm not entirely sure that this is a good idea (i.e. that can give me some significant biological insight).

@LinuxBlanket Ok, that information focused the question. Still you haven't explained what are you trying to do with this comparison. Why do you want to compare transcriptomes coming from different datasets? Having several conditions is not a problem if you want to identify what is more important of each condition. But I'm not sure what comparison are you willing to do

@llrs I'm especially interested in genes expressed exclusively in gonads, but I also want to have some insights about the basic transcriptomic features of this animal, like its housekeeping genes. That's why I need to have as much tissues as possible, so that I can really know which mRNAs are transcribed only by the gonads and which ones are transcribed by every tissue instead.

Great! (Do not forget to [edit] this into the question.) How many samples do you have from the gonads and from other tissues? If this is your objective and this is the best you can find I do not see any problem in doing it. The biological problems are hard without knowing the organism. The technical problems can't be addressed without finding them... so I'm not sure if I can answer more

As mentioned in the comments, your biological question is going to be vitally important in determining the best way to tackle this issue. However, there are some things that you should generally keep in mind.

Most de novo assemblers will perform kmer normalization. This means that the merged RNA-Seq data will be normalized across all libraries. This is almost certainly not what you want because you'll be randomly selecting reads from all libraries without regard to their treatments. So, if you're going to allow for kmer normalization, you should normalize prior to merging the libraries so that you have equal representation of the reads by library. If you don't normalize, you're going to have a very fragmented assembly, and it will take a long time to generate.
De novo assemblers using short reads assume all identified exons with junction support may be used to generate splice variants with equal probability. I'm not entirely sure how to best phrase that, but here's an example. Suppose in a cold-treated library that Gene A generates a transcript using exons 1-2-3, but a heat-treated library generates transcripts using 3-4-5 and 2-3-5. It may be that merging the libraries results in transcripts generated containing various permutations of those exons (e.g. 1-2-3-5, 2-3-4, etc.) that aren't generated in any of the biological settings but can no longer be distinguished. This will lead to false positive transcripts which will detrimentally affect expression quantification.
Merging assemblies is difficult. So, to be true to each treatment, it makes more sense to generate individual assemblies (which will also result in library-specific kmer normalization). StringTie essentially does this using reference genomes. You assemble using each library (all bio-reps for a treatment), and then merge after the fact and re-quantify. This is the best approach. However, it can be difficult and computationally taxing. Some programs like CD-HIT (specifically cd-hit-est) are typically used in these cases. Keep in mind that CD-HIT will drop sequences with high similarity, not merge them. This may lead to loss of some uniquely assembled sequence. Other approaches involve using programs like CAP3 (a DNA assembler) to treat the transcripts like longer reads for assembly. Depending on the parameter choice, this might work, but you will risk generating chimeric transcripts.

How you move forward with the assembly process will again depend on your biological questions. I would tend to promote that you tackle the questions one at a time rather than using a single transcriptome to address them all, but there may be circumstances in which merging assemblies makes more sense.
In your specific case, I would perform multiple assemblies using a variety of the approaches described above and then annotate the assemblies (ex. using Trinotate) to see how their quality metrics compare. I would discourage generating a single transcriptome and then re-assembling individual transcriptomes from that assembly -- the single transcriptome may have misassemblies that would then propagate to the individual assemblies, and this process would generally be messy/difficult.

Thank you for your answer. I edited my question to clarify it a little bit: I want to use the assembly as a reference to map reads from the single transcriptomes and do (among others) DE analysis. I probably won't use the sequence information directly from the reference.

As mentioned above, "The biological problems are hard without knowing the organism. The technical problems can't be addressed without finding them." I would second that statement and direct you to the last paragraph I just added for clarification.


text:
Retrieve ID ligand from PDB file
I have thousands of PDB files and I need to extract the following information from each of them: an ID or the name of the small molecule (ligand) present in the complex.
Is there a way to do that with BioPython or similar tools? Taking this PDB entry as an example, the tools should return 7FM.

Can you please [edit] your question and include some examples and the output you would want from them? That way, we can test our solutions using your examples and we can be sure to give you the output you need.

The PDB entry states that the ligand is 7FM.

Anyway use the MMCIF2Dict tool to extract  the information from the header via the PDB module of BioPython. If it's not in the header you can't do it because XML is not supported (not currently supported). Goodluck

And how do you use this function?

question is not clear anymore is the author referring to files a s PDB format or PDBx/mmcif format ? In case of the former Biopython can get the name of an HETATM residue see https://files.rcsb.org/view/5TO8.pdb

I found out that the package PyPDB has the function get_ligands that does exactly what I'm looking for.

Thanks, for the feedback. Apologies for not spotting this.

No problem! Strange I was not able to find this tool before...

copying from How to download ligands for PDB structure one of the answers, one of mines, code  :

import requests

import json

parsed_lig_dict = {}

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    # print(risposta.status_code)
    # print((parsed))
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    # print('\n\n',parsed_lig, type(parsed_lig))

    # for i in parsed:
    #     print('\n', i  , parsed[i])

    for lig in parsed_lig:
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")
        parsed = json.loads(risposta.content.decode())
        # print('\n',risposta.status_code, risposta.content)

        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"],
              # parsed["pdbx_entity_nonpoly"]["comp_id"])

        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] ,
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]

    for lig in parsed_lig_dict:
        for chain in parsed_lig_dict[lig][1]:
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")
            # print(risposta.status_code)

            parsed = json.loads(risposta.content.decode())
            # print('\n',risposta.status_code, risposta.content)
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])

            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]

    # print(parsed_lig_dict)

    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])

    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])

    cnt = 1
    for i in parsed_lig_dict:
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]
            chain = y
            comp_id = parsed_lig_dict[i][0]
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/{pdb_file}/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,
                              allow_redirects=True)
            print(risposta.status_code)

            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"
            open(filename, 'wb').write(risposta.content)
            cnt += 1

if __name__ == '__main__':

   pdb_file = '5to8'

   main(pdb_file)

output file : lig_7FM_B_701_1.sdf :
7FM
  ModelServer 0.9.10

 39 43  0  0  1  0  0  0  0  0  0
   -2.0610  -10.3590   40.5700 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.1170  -11.7950   42.6720 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8640  -12.6100   43.6800 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6970  -10.3290   36.5750 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8210  -10.9080   35.9890 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6600   -8.1940   35.9530 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.0300   -9.0340   41.4650 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.9390   -9.6680   41.0120 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7980  -10.0800   39.7250 N   0  0  0  0  0  0  0  0  0  0  0  0
    1.7720   -9.8430   38.8300 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9200   -9.1930   39.2330 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9910   -8.8020   40.5580 C   0  0  0  0  0  0  0  0  0  0  0  0
    4.0340   -8.8690   38.3150 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0830   -9.8810   41.9510 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3380  -10.4670   41.7760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.3020  -10.9670   40.4240 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8290  -11.6870   41.4760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8800  -11.1950   42.8370 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.0400  -12.4130   41.5840 N   0  0  0  0  0  0  0  0  0  0  0  0
   -5.2460  -12.6990   43.0180 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.5160  -10.2930   37.5490 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.4060  -11.1430   37.2090 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.8420  -10.2200   35.4050 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.7360   -8.8690   35.4000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6240   -8.9330   36.5380 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.9670  -12.3410   35.9950 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5680  -13.0870   37.1550 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.4690  -13.1930   34.5670 S   0  0  0  0  0  0  0  0  0  0  0  0
   -4.0440  -12.7310   37.3480 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.9670  -12.8650   40.6430 C   0  0  0  0  0  0  0  0  0  0  0  0
   -6.9720  -13.5050   40.9690 O   0  0  0  0  0  0  0  0  0  0  0  0
   -4.6200  -13.4130   38.5900 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.8010  -12.6200   39.1510 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5520  -12.3690   33.7940 O   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430  -14.5640   34.8900 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0530  -13.2140   33.7590 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.0490   -8.1440   38.8710 F   0  0  0  0  0  0  0  0  0  0  0  0
    4.6120   -9.9430   37.7390 F   0  0  0  0  0  0  0  0  0  0  0  0
    3.5350   -8.1080   37.3100 F   0  0  0  0  0  0  0  0  0  0  0  0
  1 16  2  0  0  0  0
  1 15  1  0  0  0  0
  2 17  2  0  0  0  0
  2 18  1  0  0  0  0
  2  3  1  0  0  0  0
  3 20  1  0  0  0  0
  4  5  2  0  0  0  0
  4 25  1  0  0  0  0
  4 22  1  0  0  0  0
  5 23  1  0  0  0  0
  5 26  1  0  0  0  0
  6 24  1  0  0  0  0
  6 25  2  0  0  0  0
  7 12  2  0  0  0  0
  7  8  1  0  0  0  0
  8  9  2  0  0  0  0
  8 14  1  0  0  0  0
  9 10  1  0  0  0  0
 10 21  1  0  0  0  0
 10 11  2  0  0  0  0
 11 13  1  0  0  0  0
 11 12  1  0  0  0  0
 13 39  1  0  0  0  0
 13 38  1  0  0  0  0
 13 37  1  0  0  0  0
 14 15  1  0  0  0  0
 15 18  2  0  0  0  0
 16 17  1  0  0  0  0
 17 19  1  0  0  0  0
 19 30  1  0  0  0  0
 19 20  1  0  0  0  0
 21 22  1  0  0  0  0
 23 24  2  0  0  0  0
 26 28  1  0  0  0  0
 26 27  1  0  0  0  0
 27 29  1  0  0  0  0
 28 36  1  0  0  0  0
 28 34  2  0  0  0  0
 28 35  2  0  0  0  0
 29 32  1  0  0  0  0
 30 33  1  0  0  0  0
 30 31  2  0  0  0  0
 32 33  1  0  0  0  0
M  END
> <model_server_result.job_id>
RFDalM1AmREvgkv0PqdfQw

> <model_server_result.datetime_utc>
2023-10-27 22:13:14

> <model_server_result.server_version>
0.9.10

> <model_server_result.query_name>
ligand

> <model_server_result.source_id>
pdb-bcif

> <model_server_result.entry_id>
5to8

> <model_server_params.name>
atom_site

> <model_server_params.value>
{"label_asym_id":"B","auth_seq_id":701}

> <model_server_stats.io_time_ms>
5

> <model_server_stats.parse_time_ms>
76

> <model_server_stats.create_model_time_ms>
2

> <model_server_stats.query_time_ms>
244

> <model_server_stats.encode_time_ms>
0

> <model_server_stats.element_count>
39

$$$$

picture of .sdf file :

my attempt, will get only non water (H2O) or non aa (aminoacid) present in the PDB .... tricky if the ligand is an aa :
from Bio.PDB import PDBParser

from Bio.PDB import Selection

d = {'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'}

pdb_file = '5to8'

struct = PDBParser(QUIET = True).get_structure(pdb_file, pdb_file+'.pdb')

res_list = Selection.unfold_entities(struct , 'R')

for i in res_list :

    if i.resname != 'HOH' and i.resname not in d :

        print(i, i.resname)

output :
<Residue 7FM het=H_7FM resseq=701 icode= > 7FM