text:
Convert SRA to FastA
I'm trying to get the FastA files for some accessions (like NC_001416.1). I did not managed to find an FTP server or direct link to these files (I want to get it from command line with wget, not from a web browser). But I found an "equivalent" file in ftp://ftp.ncbi.nlm.nih.gov/sra/refseq/NC_001416.1
This file is an SRA file that should be processed with the SRA-toolkit. In particular, fastq-dump seems to be the tools of choice.
I ran
fastq-dump --fasta NC_001416.1

and I got a NC_001416.1.fasta file. But the file contains the sequence I'm looking for in chunks of 5000 bases. The IDs in the output file are:
>NC_001416.1.1 length=5000
>NC_001416.1.2 length=5000
>NC_001416.1.3 length=5000
>NC_001416.1.4 length=5000
>NC_001416.1.5 length=5000
>NC_001416.1.6 length=5000
>NC_001416.1.7 length=5000
>NC_001416.1.8 length=5000
>NC_001416.1.9 length=5000
>NC_001416.1.10 length=3502

I need that sequence under a single ID:
>NC_001416.1 length=48502

I cannot find the proper parameters to supply to fastq-dump to achieve the desired result.
Does anyone know how could I get the FastA file I'm expecting? How to process properly the SRA file or where I could download the FastA file directly.

Out of curiosity, why do you write "FastA"? I've never seen it written that way before, only "FASTA" or "fasta". Why the final capitalized A?

No reason, really. It's a personal way of life XD FASTA comes from Fast-All, because it originally worked with all alphabets. There was also FAST-P for proteins and FAST-N for nucleotides. So FASTA is a shorthand for Fast-All. And I just used CamelCase to name it. FastA. Same for FastQ. First full word capitalized, second word reduced to a single letter and capitalized too.

@Poshi as much as I follow your logic, I find really bad not to use the "standard" .fasta. using your own capitalization generates confusion.

In file names I use the standard .fasta, .fa or .fas. In natural language, I use what it feels better for me. I'm sorry if it confuses you. To me, this capitalization is clearer. If we have to go purists, the proper way to refer the format is FASTA, all capitals (check https://en.wikipedia.org/wiki/FASTA_format or https://zhanglab.ccmb.med.umich.edu/FASTA/). But that notation hurt my eyes. It's like it's shouting to me.

typing it as FastA is kind of fun. It is a nice reminder that the file is pronounced like "fast-eh?" instead of "fast-uh".

Finally, I found an alternative to the SRA translation: a link that works! For those of you interested in knowing how to download FastA files from NCBI using an accession number, try the following link:
https://www.ncbi.nlm.nih.gov/search/api/sequence/${accession}/?report=fasta

Using wget to download the accession used as example:
wget -O NC_001416.1.fasta "https://www.ncbi.nlm.nih.gov/search/api/sequence/NC_001416.1/?report=fasta"

This works, but is it documented anywhere, and how to get a gzip?

I didn't found that in the documentation, but playing with the web. If you need a gzip, my advise is to download first and gzip later (or do it at the same time with pipes, but the gzip is on your side either way).

OK. gzipping locally is problematic because the high cost operation is the network transfer.

I'm sorry, I didn't found any other option. Try to find information about the report parameter. Maybe there is something like gzfasta. But I cannot point you to any documentation :-?

Here is a Biopython solution via Entrez utilies, simply copy the code into a python script and execute it:
from Bio import Entrez, SeqIO

Entrez.email = 'me@email.com'

def get_fasta_from_ids(ids_list):
    handle = Entrez.efetch(db='nucleotide', id=','.join(ids_list), rettype='fasta')
    for record in SeqIO.parse(handle, 'fasta'):
        yield record

SeqIO.write(get_fasta_from_ids(['NC_001416.1']), 'NC_001416.1.fa', 'fasta')

You can add more ids to the input list to retrieve multiple records at once.

Thanks for you code snippet, but this is also a solution that needs extra software to be installed (BioPython and this snippet itself), while I'm looking for something that can be run out of the box in any system with no extra software if possible.

Is the Entrez.email actually doing anything? Does the Entrez package need it to be set in order to run?

@terdon Yes Entrez.email is required (I do use my actual email address although I've never received an email about this). Entrez.api_key can be also be added to speed up queries. BTW, I agree that using the E-utilities (as in your solution and mine) is the most flexible approach, whether you do it via python or something else is mostly a matter of taste

You can either do this directly as you show in your answer or, for a more sophisticated and flexible approach, use NCBI's edirect tool:
esearch -db nucleotide -query 'NC_001416.1' | efetch -format fasta > NC_001416.1.fa

That's a quite elegant solution, but it involves installing new software, and I'm trying to generate my scripts with as less extra software as possible. Althought this is a more robust solution than having a link that can change over time.

@Poshi fair enough, I'm glad you found a solution. You might want to accept your answer to indicate that it is the one you preferred.

You can add the --origfmt option, and see what comes out. Either you will have what you want, or perhaps they uploaded the sequences that way.
See: https://ncbi.github.io/sra-tools/fastq-dump.html

With that parameter I have the exact same output :-( I think I already tried everything that made sense in my head.

Did you also try to set the --fasta to 0, so that there is no line wrapping?

Yes, same effect, oly difference is that instead of a line of ID and several lines of sequence, I have on line of ID and on line of sequence. But still 10 sequences.


text:
Multiple files for case in GDC data portal
I am downloading a set of Copy Number Variation samples from the GDC data API. I noticed that some cases have multiple files (aliquots) for the same case (see for example this case --- here the set of tumor samples with masked CNV in which it would be necessary to choose one). I need to have one file for each case in the resulting data matrix.
I looked through the GDC repository but there's no indication, the same in the GDC API.
I think they are samples processed in different institution or at different times. I looked at the metadata but I cannot find a good criteria to choose one of them. Here's the metadata of an example set for a specific TCGA-LUSC case:
{'data_release': '12.0', 'data_type': 'Masked Copy Number Segment', 'updated_datetime': '2018-08-07T15:59:14.863537+00:00', 'file_name': 'ELBOW_p_TCGAb39_SNP_N_GenomeWideSNP_6_A07_598508.nocnv_grch38.seg.v2.txt', 'submitter_id': 'e65d2f5e-96b5-4128-814e-368d830723f5_dr12_nocnv', 'file_id': '9405ea1d-1241-4e87-b1cd-a8e85008964d', 'file_size': 27971, 'id': '9405ea1d-1241-4e87-b1cd-a8e85008964d', 'created_datetime': '2018-05-20T21:25:50.066318-05:00', 'md5sum': '0faad919fb8ccfaa1cbcf04990433b2c', 'data_format': 'TXT', 'acl': ['open'], 'access': 'open', 'platform': 'Affymetrix SNP 6.0', 'state': 'released', 'version': '1', 'data_category': 'Copy Number Variation', 'type': 'copy_number_segment', 'experimental_strategy': 'Genotyping Array'},
{'data_release': '12.0', 'data_type': 'Masked Copy Number Segment', 'updated_datetime': '2018-08-07T15:59:14.863537+00:00', 'file_name': 'WHIRR_p_TCGA_168_169_170_redo_N_GenomeWideSNP_6_B10_845144.nocnv_grch38.seg.v2.txt', 'submitter_id': '9fda454b-946a-4760-8f3b-c24bdad9dbfe_dr12_nocnv', 'file_id': '91ac1284-d4e7-4be3-8266-9fec63b41cf3', 'file_size': 19377, 'id': '91ac1284-d4e7-4be3-8266-9fec63b41cf3', 'created_datetime': '2018-05-20T21:50:19.566069-05:00', 'md5sum': 'c2a55b28417de57d00c8e3eda038ca6b', 'data_format': 'TXT', 'acl': ['open'], 'access': 'open', 'platform': 'Affymetrix SNP 6.0', 'state': 'released', 'version': '1', 'data_category': 'Copy Number Variation', 'type': 'copy_number_segment', 'experimental_strategy': 'Genotyping Array'}]

In particular I thought about filtering by timestamp, but they have very close timestamps which are probably the timestamp of the final upload/update, not of the sample creation.
I noticed that one sample has the word "redo" in the filename. However, I need a methodology that can be automatized through a Python script - I cannot do it by hand.
I finally thought about taking the sample with more data, i.e. the sample with highest file size. However, this can be a bad criteria if the sample with smallest file size has a higher quality due to some kind of filtering.
Have you already encountered this problem and find a solution?
Edit
By looking at the legacy TCGA barcode, the difference is due to portions distributed on one/more plates, analyzed by the same center. The two barcodes are:

TCGA-66-2791-01A-01D-0978-01
TCGA-66-2791-01A-01D-1969-01

In boldface the plate id. My hypothesis is that by chance, the biggest one has shorter segments (more resolution, more lines) than the smallest, but both should have the same information (assuming the two aliquot have the same CNV).

have you tried tcga biolinks package available through R?  Have a look on the documentation,  may it help to solve your problem. https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html

Thanks @AmmarSabirCheema, I didn't know it - it is certainly useful!

source: Stack Exchange

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

sort: votes


text:
DESeq2 complicated design - effect of replicated samples
I have RNAseq data from a relatively complicated experimental design with variables = genotype, treatment, time, and batch. I have 2 biological replicates for each genotype/condition, however unfortunately in the first iteration of the experiment one sample had poor RNA quality and is unusable. To get around this lack of a replicate, we repeated the experiment for the two biological replicates in the condition of the lost sample, and did a new library prep/sequencing from this repeated experiment.  The repeated experiment obviously lead to a batch effect, and although it isn't the main driver of variance, it is a significant factor.  I would like to model or remove the batch effect from the experiment using DESeq2 if possible, but after reading this similar post: can I model technical replicates in DESeq2? I am concerned about inflating the power of the test for those conditions that were replicated.
Here is a diagram of the experimental design to make it more easy to understand:

My "colData" for the DESeq2 object looks like this:

My previous thought was to simply collapse the genotype+time+treatment variables into a single combined factor "condition" and then model the batch effect using a design of "~batch + condition".  However this doesn't really address the concern of inflating the power of the tests due to technical replication.
Should I Just collapse the counts of these replicated conditions even though the batch affect will affect the combined count values? Or should I use a different tool than DESeq2 ?  Or is there a way for me to remove the batch effect prior to collapsing the counts for replicates?
What would  you guys recommend?

Were the libraries in the second batch simply resequenced or were they actually recreated (presumably from the same samples, though I wonder if this is cell culture so these represent different passages).

@DevonRyan The libraries in the second batch were actually recreated from the same samples --- this is indeed a cell culture experiment, so different passage number is a concern, however these repeated samples were within ~1 passage of the samples in the original experiment (they may have actually been from the same batch of cryopreserved aliquots).  So, the experiment was repeated for the two KO cell lines and then libraries were prepared using the same techniques/kits and sequenced on the same machine.

Different passage numbers are actually a good thing in this case, since it makes them not technical replicates. I would suggest proceeding with your current design matrix.

Good point.  It is likely that there were small variations in timing of samples that also would have affected the replication.  Are you aware of any biases that could arise from this type of experimental design? In the post I linked you and others mentioned that technical replicates could inflate the power of certain comparisons. I'm specifically concerned about comparing the # of DEG's between the wild-type and KO conditions for drug vs. mock treated conditions -- if power is inflated I may see more DEG's for the KO comparisons due to more replicates right?

I agree with @devon-ryan that the design matrix looks reasonable, i.e. that it's okay to treat the new samples as biological replicates of the original batch, and that the batch number should be included in that matrix (assuming DESeq2 doesn't complain about it) to filter out most of the batch effects.
DESeq2 tends to be conservative in its differential expression calculations. While it's true that more replicates in one sample will increase the number of differentially-expressed genes that are statistically significant (because the biological variation can be better modelled), I wouldn't expect this to change the overall picture.
However, whole-genome differential expression testing should not be used on its own to demonstrate biological relevance. I recommend as a first-pass sanity check that people compare any differential expression results to what is seen on a genome browser to identify any quirky things going on (e.g. expression across only a part of the gene). Ideally, results should be confirmed with targeted analyses on the same (or similar) biological samples, either at the transcript level (e.g. qPCR) or at the protein level (e.g. Western Blot). When aggregating evidence from multiple analyses, the specific way in which one of those analyses was done is less important.
Aside: It'd obviously be better [in the future] to make sure that your sample replicates are both robust to drop-out, and allow for outliers to be noticed, reducing the batch-effect problem. In practise this means a minimum of four replicates per informative condition, but because funders / researchers are well known for only choosing the minimum, I'd just state that six replicates per sample as what is needed to produce good results.

Thanks for your comments.  I agree that it is important to validate any important findings by other methods like qPCR/Western.  And yes, I learned an important lesson about experimental design here when I had a sample drop out.  Its tough though when you want to do a time-series experiment comparing drugs and genotype because the amount of samples is x10 for each additional replicate.  Either way, I agree that in the future I will go for at least 4 replicates.  Thanks!

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


text:
Using Broad Institute RNAseq pipeline
There is a Github report for Broad's pipeline on RNA-Seq. It has a docker build file and WDL scripts.
However, I don't see how any of the WDL script is being used in the docker (Dockerfile). It looks like Broad's docker runs on the standard command-line tools, but not via the WDL language.
How do the Broad bioinformaticians take advantage of the WDL scripts?

You need to follow the readme of this folder to  use the pipeline, there you will find instructions about how to align, mark duplicates, quantificate, and aggregate the results.
If you want to execute the WDL scripts you can follow this guide using Cromwell. Basically use: java -jar Cromwell.jar run myWorkflow.wdl --inputs myWorkflow_inputs.json


text:
Understanding SNP coding for association analysis
I'm working on a project about detecting SNP association with a disease. As I understand, SNP is a single variation of the nucleotide that occurs for more than 1% of the population. So, if a gene is ATTG and there's a variation ATAG then the SNP is A. However, I couldn't connect this idea with the dataset in hand. The rows in my dataset represent each patient and the columns contain SNP information. For example:
ID  exm355 exm615
1     T_T   A_C
2     T_T   C_C
3     A_T   C_C

I have no idea why the SNP columns contain 2 nucleotides(T_T, A_T, A_C, C_C). As the definition of SNP, I thought it should show only the variant nucleotide or am I misunderstanding anything? How could I interpret T_T or C_C and how could I know which nucleotide is the variation from the common ones in the population?
The source says that the case and control subjects were genotyped by the Australian Genome Research Facility using the Illumina HumanExome BeadChip-12v1_A, which exonic content consists of $>$ 250 000
markers representing diverse populations and a range of common
conditions. An individual was duplicated to test genotyping reliability and
quality. The identity by descent (IBD) matrix between all pairs of individuals
was estimated after linkage disequilibrium (LD) pruning and used for
quality control and for the mixed linear models analyses.

You don't explain how you get the dataset, but my guess is that you are looking at the diploid genotype. Hence two nucleotides for each location (diploid genome).

If you got the genotyping data from an array experiment, one of them is a major allele and the other is a minor allele. Its rare to see multi-neuclotide polymorphisms at least in human populations whereas your data contains them at every locus. I'd venture to guess its from multiple strains of mice or other non-human organism.

Welcome to the site. Could you please [edit] the question to include the format of the file or how did you obtain it? It could also be that one nucleotide is the actual nucleotide and the other is the reference (?). Knowing more about the file would help us understand better what it means

I just add some info about the source of the dataset. Hope it could shed some light on the matter

Humans are diploid. This means they have two copies of each locus (with the exception of sex chromosomes in males). Your two alleles for each individual represent the two copies of this variant locus.


text:
How to extract Genbank sequence and identify ORFs from the sequence in R?
I used
ape::read.GenBank

and then used the function
ORFik::findORFs

to find ORFs in the sequence of interest. But the former one results in a DNAbin class format and the latter one requires FASTA format.
Is there a function available to convert the sequence from DNAbin class format to FASTA format?
(I could think of directly downloading the FASTA format file and use the function.)

You should be able to use ape's write.dna() function to output fasta format sequences


text:
Downloading full search results from BlastKOALA?
I have just started using BlastKOALA KEGG which has been useful in annotating (aminoacid) sequences. This is their website: https://www.kegg.jp/blastkoala/
When you get results, there are links for downloading. However these links will not download all detailed query search results, but just general notes already on the screen. To get these results I need to manually click on each query result on the page, which becomes impracticable with >500 entries. Thus I think need is a tool to download all linked contents from a webpage. I have been trying 'wget' however it doesn't work. It says 'Requested Job Not Found' whatever I do.
Please, did anyone every try to achieve this? Thanks in advance.
P.S. I have provided a working link as an example, below:
https://www.kegg.jp/kegg-bin/blastkoala_result_list?id=7345faab3969a3ef10a4a7542ba749386068138f&passwd=NKfm46&type=blastkoala&code=user
I think it should remain active for a week.
I have inserted a picture showing examples of query result links below:

At this moment I would like to end up with the listed hit results for each of the queries, whose links are showing in blue; one result per file. In short, it is using wget for downloading all links in the page into a dedicated folder.
Thanks in advance!

Could you provide some links of the entries you want to download and the command with wget you tried? This would be useful to create a program or change the command line.

@Llopis Sure, I'm adding that now to the question. I think they keep temporary links for results. Thanks!

I can't see any link to the entries you want to retry on the link you posted.  But you can download the annotations. Perhaps is that what you want?

Yes, please [edit] your question and show us an example of the data you want to download. I don't understand what you need from your question.

@terdon thanks, I have done just that. I may have found a workaround, but if anyone provides a better solution it would be great.

So you want to download the target of each of the links in the "Query" column? Do you also need the "Second best" ones? And how do you want them organized? One file per "query"? Everything in one file? Do you just want the "gene (Length)" column from the target? Please [edit] and give us an example of what you want to end up with.

Perhaps it might be easier for you to explain why your current "workaround" is not good enough for you.

@Llopis Following scattered suggestions I've gathered from related download threads I am using the command line wget -r -l3 -H -t1 -nd -N -np -erobots=off "link_path" which is returning all hits, but in a very disorganised way and raw format. I am no programmer myself.

How would it be better organized? Not being a programmer doesn't prevent you to program (or to learn to do it :) and even less to explain yourself. Please continue clarifying the question by answering terdon answers


text:
What is the correct way to map Hi-C data with bwa mem?
Library Prep

I have a Hi-C library prepped using an enzyme that cuts at GATC, so it leaves GATCGATC as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.

Data Pre-processing

The reads were adapter and quality trimming with trimmomatic paired-end mode, keeping reads that overlap.
Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (GATCGATC). For example, the read:

5' AAAATT**GATCGATC**AATTATGGATAA 3' becomes: 5' AAAATT**GATC 3'

In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything 3' of the junction sequence was removed.

My Question

What is the correct way to now map these trimmed Hi-C reads to a genome using bwa mem?

Progress so far/things I have tried

In the past I have just tried mapping using bwa mem with no special parameters using this command: bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz. One concern is that the computations that bwa mem performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).
Another option to consider is mapping each file, hi-c_f.fq.gz and hi-c_r.fq.gz, independently then adding mate information with samtools fixmate. However, I do not know how to do this.
A last option that I have considered is to use bwa mem in paired-end mode, but to use a special flag. In the past I believe that -5 was the correct flag, but in this post Heng Li discusses that the -5 flag behavior has changed, as well as a new -q flag. In this post Heng Li discusses that he thinks bwa mem is still a good option for mapping Hi-C data but does not discuss the proper flags.

Limitations

The answer should provide the correct way to map Hi-C using bwa mem, but other software suggestions will also be considered.

This ends up depending on the software you want to use downstream of mapping. For example, HiCExplorer still needs reads to be aligned as single-end. I imagine it's not alone in that regard.

See here, in particular slide #10 from the tutorial:
bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.

Thanks for your post. It looks like the -SP flags are sufficient to map the trimmed data as I discussed above, while -5SPM is appropriate for Hi-C data that do not have the linker sequence trimmed.

There are some interesting suggestions from the HiCExplorer developers:

Mates have to be mapped individually to avoid mapper specific
  heuristics designed for standard paired-end libraries.

# map the reads, each mate individually using
# for example bwa
#
# bwa mem mapping options:
#       -A INT        score for a sequence match, which scales options -TdBOELU unless overridden [1]
#       -B INT        penalty for a mismatch [4]
#       -O INT[,INT]  gap open penalties for deletions and insertions [6,6]
#       -E INT[,INT]  gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] # this is set very high to avoid gaps
#                                  at restriction sites. Setting the gap extension penalty high, produces better results as
#                                  the sequences left and right of a restriction site are mapped independently.
#       -L INT[,INT]  penalty for 5'- and 3'-end clipping [5,5] # this is set to no penalty.

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R1.fastq.gz 2>>mate_R1.log | samtools view -Shb - > mate_R1.bam

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R2.fastq.gz 2>>mate_R2.log | samtools view -Shb - > mate_R2.bam

# build matrix from independently mated read pairs
# the restriction sequence GATC is recognized by the DpnII restriction enzyme

$ hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
                 --binSize 10000 \
                 --restrictionSequence GATC \
                 --threads 4
                 --inputBufferSize 100000
                 --outBam hic.bam \
                 -o hic_matrix.h5
                 --QCfolder ./hicQC

Option -SP exactly aims to disable "mapper specific heuristics designed for standard paired-end libraries".


text:
How can the Autocovariances, autocorrelations, and autocorrelation coefficients be calculated from a Protein Amino Acid Sequence?
Given a normal protein sequence with the 20 standard amino acids, how can the 'Autocovariances', 'autocorrelations', and 'autocorrelation coefficients' of the sequence be calculated?
What is meant by these terms in the context of a protein sequence?
What are the appropriate algorithms and formulas?
Is there any existing package/function to calculate these values?

Could you provide some context? It could clarify the question and attract more answers. Maybe you read that in an article or it is for some specific goal in mind, or someone asked you this because a discussion about another topic

Protein autocorrelation is similar to usual autocorrelation of time series (that is, function $r(k)$ measuring correlation between values lagged by $k$), except that instead of time dimension we slide along the backbone of amino acid chain. The "dependent variable" can be any quantifiable property of the amino acids - then the autocorrelation measures similarity in that particular property.
See e.g. Fernandez et al., 2008, for more extensive background:
https://www.tandfonline.com/doi/full/10.1080/08927020802301920?scroll=top&needAccess=true


text:
RNA velocity: competing explanations for variable ratio of spliced to unspliced transcript
In "RNA velocity of single cells", La Manno et al. look at ratios between spliced and unspliced mRNA as a way of estimating the velocity of a cell through transcriptome space.
In checking for a universal steady-state spliced-unspliced ratio, they write "We found that the steady-state behaviour of most genes across a wide range of cell types was consistent with a single fixed [ratio between spliced and unspliced] (Extended Data Fig. 3a–c). However, 11% of genes showed distinct slopes in different subsets of tissues (Extended Data Fig. 3d, e), suggesting tissue-specific alternative splicing (Extended Data Fig. 3f) or degradation rates."
Tissue-specific alternative splicing or degradation rates certainly could account for these patterns, but how do La Manno et al. rule out that some tissues are simply not in a steady state? The figure they cite (screenshotted below) uses data from the Tabula Muris adult mouse cell encyclopedia. In an adult mouse, not all cells would not be expected to be in steady state -- for example, T-cell progenitors from the thymus and epithelial cells in many organs undergo continual replenishment.

If cells aren't in a steady state then one wouldn't expect to see the results in (e), where there are two slopes. One would instead expect to see results more like Plekha3 in panel (c), where there's more variability around the the steady state line (or perhaps no steady state line at all).

Could you explain why that is that case? For example, couldn't you get two slopes as follows? Suppose that keratin 17 increases during development of multiple epithelial lineages, but is constant in non-epithelial cell types. For epithelia, there would be an excess of unspliced transcripts due to the portion of the population that is developing. The epithelia would have a line with a high slope and the rest would have a line with a low slope.

The difference between developing and differentiated is a gradient.

I think I understand. So if cells are on average less differentiated in intestinal epithelia and more differentiated in tracheal epithelia, you'd expect them not to fall a line of the same slope?

You wouldn't expect undifferentiated cells to fall on any particular line, since that's a broad spectrum of states. They may have some line that they tend to fall on, but you'd then expect to see a lot of noise between that and the "differentiated line", due to various levels of differentiated status.

In addition to Devon Ryan's data-driven answer, I've found a clue in the text. La Manno et al have not used the full Tabula Muris: the figure legend indicates that they used only terminally differentiated cells.


text:
GATK documentation for required depth to reliably call heterozygous mutation in diploid organism?
I'm looking for official GATK documentation (or a recent manuscript) that defines a general recommendation/requirement for sequencing depth to reliably call a heterozygous point mutation in a diploid organism (WGS). In this case, I'm working with humans. Specifically, suppose we're looking at position X in the genome and we want to classify (True/False) whether there was enough coverage to detect a heterozygous point mutation.
What coverage would that be? I'm curious about both: (1) UnifiedGenotyper and (2) HaplotypeCaller
I know general practice is to have 30x coverage genomewide, but I don't see anything in GATK's documentation (i.e., Best Practices) about required coverage.
Really appreciate your help.
Note: I also submitted a question to the GATK team end of last week.

UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement:

We believe the sites around 30X (and above 150X) are being filtered
  due to low mapping quality (since it is whole genome so there are many
  areas that are hard to map), which explains the low sensitivities in
  these areas. As the error bars show, there are very few variants in
  these low mapping quality areas.

She then posted the following plot that falls much more inline with what we would all expect. This is a plot showing sensitivity for NA12878 compared to NIST truth.

shlee, from the GATK forums, pointed me to a poster from AGBT 2018 where the authors compared the HiSeqX and NovaSeq instruments. Figure 3 shows sensitivity for SNPs and short INDELs across multiple allelic fractions (AF) from spike-ins. I was surprised to find that SNPs where AF = 0.4 only had ~50% sensitivity at ~25x coverage, and only ~75% sensitivity at ~50x coverage.
The poster is targeted at somatic variation, but data from figure 3 should behave like a standard diploid germline variant because you often see a wide range of allelic fractions from germline variants. They used spike-ins for figure 3. It's not clear to me whether they used HaplotypeCaller in figure 3, but I don't see how they could have used Mutect2, like they did in figures 1 & 2.
Based on this information, it appears the standard 30x genome-wide coverage is generally insufficient to confidently call many heterozygous mutations, and that ~90-100x is ideal.
Am I missing something? Are we really missing that many heterozygous mutations in standard WGS/WES?
I'll leave this question unanswered for now to allow others to chime in on interpretation, etc. I'd appreciate additional thoughts.


text:
Seurat with normalized count matrix?
I know that in Seurat we have the function CreateSeuratObject from which the analysis starts, but it accepts raw count matrix according to the documentation. I have only the already normalized count matrix, so is there a way to work with Seurat using normalized data? I would try maybe to feed in the normalized data instead of raw to CreateSeuratObject, but I do not know what Seurat is doing behind the scenes, whether it is a good idea. Besides, Seurat provides by default only one log - normalization method, but I may want to normalize the data by myself with various methods and only then start the analysis with Seurat - that is the other reason why I want to find a way to start from normalized data.

This was addressed by the Seurat developers here:

if you have TPM counts, I suggest you don't use
  Seurat::NormalizeData(), since TPM counts are already normalized for
  sequencing depth and transcript/gene length. Note that
  Seurat::NormalizeData() normalizes the data for sequencing depth, and
  then transforms it to log space. If you have TPM data, you can simply
  manually log transform the gene expression matrix in the object@data
  slot before scaling the data.

And here:

If you want to run your own normalization function, please place the
  data in log scale before placing it into object@data

And here:

If you're supplying already normalized values to CreateSeuratObject,
  you should skip the normalization step. You will likely see a warning
  when you run ScaleData (which is checking whether you normalized using
  Seurat). You can ignore this or turn it off by setting check.for.norm
  to FALSE.

This is a lot of quotes, but I found them all to be helpful in different ways.


text:
How do population genetics people define a population?
How do population genetics people define a population?
Do they define it as a layman will do? say Africans, Americans, or so? Or is there a more scientific way of doing so? For example, I think defining one's population as one's allele frequencies makes more sense, but has this point been discussed previously?

Do you talk about a general definition or you want to what what defines human populations?

It really depends on the study and the ethnicity. Mostly they use self-defined ethnicity and often they will require that three or four grandparents are of that ethnicity. Here's some detailed info on how the ethnicity was defined for the 1000 Genomes populations. You'll see that some populations (such as ASW) require all four grandparents to identify as that ethnicity, but some (such as ACB) only need three, and the cultural sensitivities in Japan meant that JPT participants weren't questioned directly but told that the study was only interested in people with Japanese grandparents. The mixed American super-population (AMR) is an exception to this rule due to the extreme admixture of ethnicity in these regions, and is based more in geography than ethnicity – this means that the individuals had to be from the place and have grandparents from the place.

A population is defined as a group of people (or organisms of the same species) living together in a geographical region. However, for human genetic studies, it is wise to use ethnicity as well (so group of people with the same ethnic background). See here how they define population in the 1000 genomes project.


text:
TagReadWithGene missing when using latest version of Drop-seq_tools
I am using Drop-seq_tools to analyze scRNA-seq in a similar way of this paper. I am using the same data (GSE97930) to get used to Drop-seq_tools.
On the GEO repo the authors provide the relevant code (GSE97930_README.txt). However, one tool used in the script is missing: TagReadWithGene.
The latest version of Drop-seq_tools does not have such method (as it should be provided by this library), and there is no version history for the whole package, so I cannot download previous versions.
The best matches for the missing method are TagReadWithGeneExon and TagReadWithInterval.
In the methods of the cited paper they say:

Reads that mapped to intronic or exonic regions of genes as per the GENCODE gene annotation were recorded.

I think that in this case the tool to use should be TagReadWithInterval since the other one should consider only exons.
I cannot find any other documentation about the tools. Do you know/already encountered this problem?

As per this response from Drop-seq_tools team, at the moment intronic regions cannot be annotated and TagReadWithGeneExon should be used instead.
However, in the version of the public repository:

TagReadWithGeneExonFunction is the old code, for people who want the more simple tagging methodology and rely on it for 3rd party tools, and TagReadWithGeneFunction, which for your purposes replaces TagReadWithGeneExon. [...]  You’ll see additional options to extract expression for CODING, UTR, INTRONIC reads [...]


text:
Histosketch vs Count-Min sketch: preserving similarity
A preprint describing a new tool and its application to microbiome analysis was recently published in bioRxiv[1]. At the core of this new tool, HULK, is a new data structure called a histosketch[2] which is similar in spirit to CountMin sketches and MinHash sketches that have become quite popular in bioinformatics in the last several years.
One of the main advantages the histosketch authors claim over CountMin sketch is that CountMin sketches do not preserve similarity. From section II of their paper[2] (emphasis mine):

It is practically difficult to compute histograms for data streams with typically unknown cardinality and which thus require and unbound amount of memory to maintain the histogram. In this context, the count sketch and the count-min sketch and other online histogram building methods were proposed to approximate the frequency table of elements (i.e., histograms) from a data stream with a fixed-size data structure. However, the resulting sketches do not preserve similarity between different data streams.

In what sense does the CountMin sketch not preserve similarity between streams? Does this refer to similarity at intermediate states (while data is still streaming) or does it refer to inherent characteristics of the data structure itself?

[1] Rowe WPM, Carrieri AP, Alcon-Giner C, Caim S, Shaw A, Sim K, Kroll JS, Hall L, Pyzer-Knapp EO, Winn MD (2018) Streaming histogram sketching for rapid microbiome analytics. bioRxiv, 408070, doi:10.1101/408070.
[2] Yang D, Li B, Rettig L, Cudré-Mauroux P (2017) HistoSketch: Fast similarity-preserving sketching of streaming histograms with concept drift. Proc. IEEE Int. Conf. Data Mining, doi:10.1109/ICDM.2017.64.

Answer from @user172818, converted from comment:
My understanding is that you can't compute Jaccard similarity given two count-min sketches. You can with two histosketches, though I don't know exactly how yet.


text:
Check if a set of genes fall more on a specific cluster than other clusters
I have a network which I clustered into 11 smaller clusters. I have a set of genes and I would like to check if statistically significant number of these genes fall on any specific cluster. One method which I thought of was to find which percentage of a total cluster were made up of my genes of interest. For example, say 60% of the genes in a cluster were the genes in my list. However, I was hoping that there was a better way to do this. This study on honey bees has done something of the kind.
I would appreciate if there was a way to do something similar. I do not know how to replicate what they have done in the study. Any help is appreciated.

Do you have problem with the hypergeometric test or with the later methods they do?

@Llopis I have never done hypergenometric tests before and I am not sure how to do it.

A hypergeometric test is related to the fisher test, in a table with one  dimension being some factors and the other one some other factors it is tested if they are independent. In R you con do this:
data <- matrix(c(0, 1, 2, 5), ncol = 2, nro = 2)
fisher.test(data)

Or you can use the specific function phyper, to see the relation between them you can read this question and answer


text:
What is and how to detect a dephased read
In the methods of this paper, the authors say:

To simplify analysis, we first removed any dephased reads in our library (last 6 bases of read did not match the expected sequence).

I have read this post and this one, and I think that a dephased read is one (according to the first paper) where the last 6 bases are not detected with good quality. In particular, I think this may be an example of a dephased read:
@...
CACACACACCAGTAAAAAANAGTTGCNNCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
AAAAAEEAEEEEEEEEEE6#EAEEEE##E#####################################

where the quality of the read drops in the middle. Is that correct? If yes, should they be all of bad quality or at least one? And, if not, how can I detect a dephased read?

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be.  They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases.  This would include N's, but a small deletion or insertion would also cause the last 6 bases to not match.  This is in theory detectable and fixable; the dropSeq protocol for instance catches instances where improperly short cell barcodes shift the UMIs, but these authors look like they decided not to try and bother, and just throw out everything that didn't match their exact expectations.

Interesting, thanks. Is the dephased read considered only in the reverse read (R2), where barcodes and UMI are sequenced? I assume that during the mapping of the forward read the aligner is able to deal with dephased reads - while during barcode extraction dephased reads may introduce issues, if I understood correctly your answer.


text:
Source for whole genome comparisons at NCBI Genomes
The NCBI Genomes database has these dendrograms for (presumably) whole genome comparisons for certain species, e.g. Pseudonomas aeruginosa  or Escherichia coli.
How were these comparisons done? Someone knows the source / paper?

They say "based on genomic BLAST". No useful information can be found, though.

I inquired into the details of the dendrograms after becoming frustrated with the lack of information. As with Ensembl, I'm sure that the folks at NCBI have a standardized pipeline that they run the sequences through to generate these dendrograms as no specific source deserving attribution was mentioned:

The tree is based on a pairwise, BLAST comparison of chromosome sequences (an assembly is used for genomes with 2 or more chromosomes). BLAST identity is normalized on sequence length. The distance map is generated by a Neighbor Joining algorithm, where the distance is the blast score.
Dr. Wayne Matten, NCBI

Thanks for inquiring. Too bad that their response doesn't contain more details or a link to the pipeline. Seems like they calculate the average PID of all alignments reported by blastn und use that as a distance measure. I wonder if they also consider the parts of each genome pair that are not covered by alignments.


text:
How to detect mismatch before mapping in RNA-Seq data
In the methods of this paper, the authors say:

Reads were then filtered based on quality score in the UMI region. Any read with >1 low-quality base (phred <=10) were discarded. Reads with more than one mismatch in any of the three 8 nt cell barcodes were also discarded. The cDNA reads (Read 1) were then mapped to either a combined mm10- hg19 genome or the mm10 genome using STAR (44).

I understand how to filter based on mismatch after mapping but in this case it is done before mapping. Has this to be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded"?

At each step in the SLiTing and pooling the barcodes added are not random, but are taken  from a predefined list of possible barcodes.
When the reads arrive, first you filter by reads that have >1 low quality bases in the UMI region (which is truely random).
Then you compare the sequences in the 3 cell-barcode regions against the list of allowed barcodes for each of those positions. If a read is not within an edit distance of one from one of those barcodes, it is discarded.
Finally the reads that pass these filters can be aligned.

It should not be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded". They only care if the cell barcodes match what are valid cell barcodes within 1 mismatch, the individual qualities are ignored. This sort of filtering turns out to be pretty common in practice, since you have to write a maximum edit distance into whatever you use to demultiplex by cells (or labeling each read, as was likely done in this case).

Cell barcodes =/= sequence reads.  The default settings on Illumina's bcl2fastq work exactly the same way; they allow a read to count for a sample if there is a single mismatch in the index, but not more.


text:
The visualisation of a list of genes on URD object
Seurat R package has some functions like FeaturePlot, DimPlot and DoHeatmap by which we can plot the expression of a list of genes on cell clusters. I am working with URD that likely does not have such options or I can not find that. URD has plotDot but would be clear only for a few number of genes (please look at the picture). Could somebody please help me in plotting the expression pattern of a list of genes in URD? As I am not good in R at all to write my own function or adopt Seurat functions to URD

For example Seurat DoHeatmap gives such a picture by which we can see these set of genes are being expressed on which group of cells
DoHeatmap(object ,genes.use=a vector of genes)

Or we can plot a set of genes on seurat feature lot to they are being expressed more on which cluster
FeaturePlot(object = pbmc_small, features.plot = 'PC1')

In URD only by
plotDot(object.6s.mnn, genes = a vector of genes, clustering="Infomap-60")

We can plot a set of genes but as you are being in the first picture that would not be clear much
My question is how I can plot a vector of genes on URD object
for example Seurat DimPlot function can visualise the expression pattern of a group of genes on cell clusters but when I am trying to do the same with URD, returns this error
> plotDim(axial.6somite, c("DDB_G0291149","DDB_G0291029"), plot.title="SHHA (Notochord marker)")
Error in names(data) <- cells.use :
  invalid to use names()<- on an S4 object of class 'dgCMatrix'
In addition: Warning messages:
1: In if (label.type == "meta" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
2: In if (label.type == "group" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
3: In if (label.type == "sig" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
4: In if (label.type == "gene" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
>
>

Actually only one gene can be visualised by this function on URD

Which plot do you want? Could you post an example from the similar function from Seurat? You say you are not good at R, did you try learning with a bootcamp or with online books (such as this one?

Please [edit] your question and give us a minimal, complete, verifiable example that we can use on our machines to reproduce your code. Those users who can answer will need to have access to a test dataset they can use to test their answers. Not your entire data, just a small part of it, or a similar dataset that we can play with.

Sorry, I know something in R but not as good as to be able write a function in a deadline

I just added URD clustering labels for each cel to seurat object and then I did DoHeatmap by using group.by="URD clustering labels"


text:
How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`?
Very recently a BLAST parameter -max_target_seqs n got a lot of attention. Instead of the intuitive interpretation (return the best n sequences) the parameters asks blast to return the first n sequences that pass the e-value threshold. This is affecting thousands of workflows and analyses that assumed the intuitive interpretation.
I also found that there is another parameter called -num_alignments. What is the difference? Related question on BioStars seems that is misinterpreting the max_target_seqs parameter. Is num_alignments also misunderstood?
Now the main question. How to actually run blast to give me the best n hits in the database?
Till now I thought that this example will do the job:
blastp -query $QUERY -db $PROTEINS -out $BLASTOUT -evalue 1e-10 -outfmt 6 -num_alignments 5

-- edit --
Now I understood that i can not achieve the desired output only by parameters of blast, but I still don't see an answer to how can I get the best n hits. So I wanted to check out how to get all the results and that's not intuitive either. The parameter -max_target_seqs is set by default to 500. Does it mean that if there are more than 500 significant hits, I don't have a guarantee for the best one? Do I have to specify max_target_seqs to some crazy high number to be sure that I got them all?

relevant https://bioinformatics.stackexchange.com/questions/2846/why-does-a-very-strong-blast-hit-get-lost-when-i-change-num-alignments-num-desc

To your last edit: I would think so, without knowing the source code it might be needed. However, I think that the default value is the same as what they use for the web, so if you want to have the same results on your inline command and on the web it might better to leave it as is. (In case you wanted to modify the c++ code I tried to see the license of the software, but I couldn't find anything)

The -max_target_seqs is also related to -evalue which does not work as one would think. Looking at the blast news one finds that since version BLAST 2.2.27+ released, one should use the -max_target_seqs:

4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max_target_seqs to control the output rather than –num_descriptions and –num_alignments.

So it seems that num_alignments is also misunderstood. I came to the conclusion that one cannot expect to get the best n alignments, unless you get all and filter yourself the best ones.

For a manual filtering approach look at my answer.

As a follow up to @Llopis' answer (no native blast solution), I wrote a script blast_filter.py and I added it into my collection of useful scripts generic_genomics. The script filters out the blast output and it works with streams. To get the n (5) best blast hits by e-value you can run:
blastp -query $QUERY -db $PROTEINS -evalue 1e-10 -outfmt 6 | blast_filter.py 5 > $BLASTOUT

I assume that the blast output is sorted primarily by query name and secondarily by e-value. Also, note that default max_target_seqs in blast is 500, if there will be a query with 500 alignments, script will drop a warning (the query has no guarantee of the best hit found).


text:
Reference genome for allele specific expression
We are trying to sort out a pipeline for doing allele specific expression. Our plan is to call SNPs from RNA-seq data and combine with known SNP annotations. A well known problem in ASE is reference bias, where reads are more like to map if they have the reference sequence than the alt sequence.

We have come accross two possible ways to deal with this. In many papers (e.g. Zhuo et al, Sigurdsson et al) a new reference is generated, substitute N bases where there are SNPs using FastaAlternativeReferenceMaker from GATK and reads are remapped to that.
HISAT2 offers the possibility to provide SNP calls to the index builder, and build an index that accounts for these SNPs.

Does anyone have any advice on which of these might be better?

As you said, this is a well-known problem for ASE. In Degner et al. 2009, they state that, "Perhaps surprisingly, we found that masking known SNPs does little to eliminate inherent biases in read-mapping."
As such, substituting N bases is not likely to completely solve your issue. Depending on the aligner, these N bases may also negatively affect your ability to map reads.
Of your two suggestions, I believe HISAT2 will be the better approach. There are also more intensive methods you can use such as empirical Bayesian methods or more traditional variant phasing (ex. AlleleSeq).


text:
Block wise protein imputation
I am currently working on a dataset that contains 50 samples (10 samples * 5 blocks). The features of the date set are:

The data is perfectly balanced between blocks, with equal treatment representation in each block. Each block contains 2 control samples (CTL) that are pools of all other treatment samples. These are technical replicates that are identical across multiplexes.

There is a lot of missing data (~40%),

most of the missing data is missing by block, meaning that the protein is present in some blocks, but not others.

When looking at the expression of these missing proteins in the other blocks, they tend to be around the median value for the experiment, which suggests to me that these values are missing completely at random.

I could proceed using complete cases, but then I'm missing out on a large portion of data. I am considering using an imputation method, but I'm concerned that (due to the nature of the missing data) that I will just be introducing a block effect into the data.
Any input on potential methods for imputation would be greatly appreciated.
Experimental design:
SampleNumber    Multiplex   Treatment
1   A   CTL
2   A   CTL
3   A   TMT1
4   A   TMT1
5   A   TMT2
6   A   TMT2
7   A   TMT3
8   A   TMT3
9   A   TMT4
10  A   TMT4
11  B   CTL
12  B   CTL
13  B   TMT1
14  B   TMT1
15  B   TMT2
16  B   TMT2
17  B   TMT3
18  B   TMT3
19  B   TMT4
20  B   TMT4
21  C   CTL
22  C   CTL
23  C   TMT1
24  C   TMT1
25  C   TMT2
26  C   TMT2
27  C   TMT3
28  C   TMT3
29  C   TMT4
30  C   TMT4
31  D   CTL
32  D   CTL
33  D   TMT1
34  D   TMT1
35  D   TMT2
36  D   TMT2
37  D   TMT3
38  D   TMT3
39  D   TMT4
40  D   TMT4
41  E   CTL
42  E   CTL
43  E   TMT1
44  E   TMT1
45  E   TMT2
46  E   TMT2
47  E   TMT3
48  E   TMT3
49  E   TMT4
50  E   TMT4

The data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression.
The controls are technical replicates, but they are pools of all other 40 samples. This is for TMT proteomics, where the pools are used by the MS to identify peptides across the multiplexes, not necessarily for multiplex correction. I've included multiplex (as well as treatment) in the linear model to identify changes associated with the treatment, while removing the effect of the multiplex. I have confirmed (with PCA and PVCA) that the multiplex effect is negligible after including it in the lm. My problem at the moment is the missing data.
I calculated the missing percentage from the data matrix by taking the number of missing values divided by the total number of potential values. he data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression. To calculate the missing values, I divided the total number of NA's by the total number of values in the matrix (number of samples * number of proteins). Missing data is common in proteomics, but normally this can be attributed to something, most often being close to the limit of detection. But that is not the case here, if a value is missing, it is missing in the whole multiplex, including the pool, where we know it should be found.
The missing values in this case are random, and this is definitely due to the MS. It will select the first ~65k peptides that it identifies to trigger an additional MS event, so if it doesn't see that first peptide, it will be missing from the entire multiplex. There are methods for imputation of data that is missing at random, but I'm concerned that it will just be imputing by multiplex.

40% missing data is huge. Missing data analysis is complicated on the underlying distribution. If the data set is periodic then missing data periodicity is needed. Non-periodic data can be solved using regression, I think SciPy has an automated method for this, but the volume of missing data makes this approach complex. If you can bring the missing data down to 20% by compartmenatlising the data set that would be the starting point.

This looks like a TMT proteomics experiment. What tools/pipeline did you use to arrive at your quantitative values?
There are a variety of tools that are useful for analyzing this type of data and supply appropriate imputation methods. If you are comfortable programming in R, one I would suggest is the MSstatsTMT R package: http://msstats.org/msstatstmt/
That being said, any ways you can improve the upstream analysis to decrease the amount of missing values will be extremely beneficial.


text:
Calculating bit score: How do you find lambda and K?
To calculate bitscore from score you can use this equation: $S' = (lambda*S - ln(K)) / ln(2)$
If I am trying to manually calculate the bitscore of an HSP of a pairwise blastn alignment, and I know the alignment score, how do I calculate the parameters K and lambda?
Is there a more simple equation for pairwise nucleotide alignments?

NCBI BLAST uses tables of pre-computed values for lambda/K for gapped protein and nucleotide comparison and their values depend on the amino acid substitution matrix or the substitution scores + gap penalties.
You can find these tables (from line 589 for blastn) in the file c++/src/algo/blast/core/blast_stat.c from the blast source package.


text:
Circos plot issue
I have taken out the data like chromosome number, start, and end coordinate then gene name but the circos need the data in a specific format as given in their tutorial circos
So far I have my data I have collected is as such
Symbol  Chr Start   End
RBM11   hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21 14216130;14216145;14216153;14216154;14216178;14219553;14219553;14219563;14219563;14221097;14221097;14221097;14221097;14221097;14224375;14224438;14224438;14224438;14224438;14226859;14226880;14226880;14226880;14226880 14216282;14216282;14216282;14216282;14216282;14219725;14219725;14219725;14219725;14221453;14221169;14221169;14221169;14221169;14224537;14224537;14224537;14224537;14224537;14227054;14228372;14227384;14228173;14228372
ELAVL4  hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1 50048014;50103910;50105690;50106292;50108933;50108966;50109626;50109706;50109816;50141772;50144953;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177130;50179544;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50200809;50200809;50200851;50200851;50200851;50200851;50200851    50048182;50104000;50105950;50106393;50109198;50109198;50109964;50109964;50109964;50141858;50145197;50145178;50145197;50145197;50145197;50145197;50145197;50145184;50145197;50145197;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50181300;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50197585;50197467;50197467;50197467;50197467;50197467;50197467;50197467;50203054;50201402;50202294;50203786;50201399;50203770;50201868
PTRF    hs17;hs17   42402452;42422627   42405388;42423517

So i have to make the the data as such
chr start end value [options]
...
hs3 196000000 197999999 71.0000
hs3 198000000 199999999 57.0000
hs4 0 1999999 28.0000
hs4 2000000 3999999 40.0000
hs4 4000000 5999999 59.0000

As of now i have multiple
RBM11   hs21;hs21;hs21;hs21;hs21;hs21

in for a single gene so how do i make it one ?or
RBM11   hs21

Any help or suggestion would be highly appreciated
I tried:
 sed -r ':a; s/\b([[:alnum:]]+)\b(.*)\b\1\b/\1\2/g; ta; s/(, )+/, /g; s/, *$//'

Updated answer so far
  library(splitstackshape)
df <- LEUKEMIA_CIRCOS_gene
Chr <-unique(cSplit(df, "Chr", sep = ';', direction = "long"))

this gives me this output
RBM11
hs21
14216130;14216145;14216153;14216154;14216178;14219553;14219553;14219563;14219563;14221097;14221097;14221097;14221097;14221097;14224375;14224438;14224438;14224438;14224438;14226859;14226880;14226880;14226880;14226880
14216282;14216282;14216282;14216282;14216282;14219725;14219725;14219725;14219725;14221453;14221169;14221169;14221169;14221169;14224537;14224537;14224537;14224537;14224537;14227054;14228372;14227384;14228173;14228372
2
ELAVL4
hs1
50048014;50103910;50105690;50106292;50108933;50108966;50109626;50109706;50109816;50141772;50144953;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177130;50179544;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50200809;50200809;50200851;50200851;50200851;50200851;50200851
50048182;50104000;50105950;50106393;50109198;50109198;50109964;50109964;50109964;50141858;50145197;50145178;50145197;50145197;50145197;50145197;50145197;50145184;50145197;50145197;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50181300;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50197585;50197467;50197467;50197467;50197467;50197467;50197467;50197467;50203054;50201402;50202294;50203786;50201399;50203770;50201868

my chromosome name are sorted out but the start and end remains because if I make concatenate the column I want to separate I get a really long list
    Start <- unique(cSplit(df,c("Chr","Start","End"),sep = ';', direction = "long"))

RBM11
hs21
14216130
14216282
2
RBM11
hs21
14216145
14216282
3
RBM11
hs21
14216153
14216282
4
RBM11
hs21
14216154
14216282
5
RBM11
hs21
14216178
14216282
6
RBM11
hs21
14219553
14219725
7
RBM11
hs21
14219563
14219725
8
RBM11
hs21
14221097
14221453

Thanks for providing a nice set of data to test with. However, could you tell what have you tried ans why it didn't satisfy it? We can then avoid suggesting the same route or improve your existing path.

I tried something like this from a solution in stack  sed -r ':a; s/\b([[:alnum:]]+)\b(.*)\b\1\b/\1\2/g; ta; s/(, )+/, /g; s/, *$//'  but i see trails of hs20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;, after the hs20 as an example how to resolve that?

Well you seem to try to use a single regexp to solve a problem that needs to divide a line into several ones. I would count the ; and then pair each value of each column into a new line.

okay let me try to get around if i can do it myself ..

That's the spirit! Don't hesitate to come back if you have problems, or meanwhile perhaps I manage to answer it.

partially managed to do it through  a R library..i will update my question and you can see and suggest..

How are you reading this file? It seems like tab separated, but you seem to relay on cSplit to split by columns


text:
Spearman correlation between two genes
I have FPKM data for three genes like below:
     Samples      gene1 gene2 gene3
TCGA-2Y-A9GS-01A   0.9   7.4  1.0
TCGA-2Y-A9GT-01A   0.8   1.0  0.3
TCGA-2Y-A9GU-01A   0.6   2.0  0.2
TCGA-2Y-A9GV-01A   1.2   0.5  0.1
TCGA-2Y-A9GW-01A   3.8   2.1  0.4
TCGA-2Y-A9GX-01A   2.3   2.0  1.5

I used cor.test
cor.test(~ gene1 + gene2, data = df2, method="spearman", continuity=FALSE, conf.level=0.95)

    Spearman's rank correlation rho

data:  gene1 and gene2
S = 5686100, p-value = 6.083e-09
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho
0.2984045

I have a warning message which I didn't see before.
Warning message:
In cor.test.default(x = c(0.9, 0.8, 0.6, 1.2, 3.8, 2.3, 3.8, 0.4,  :
  Cannot compute exact p-value with ties

Do I need to care about this warning message? Is it good using FPKM data for correlation?
For plotting I used ggscatter.
ggscatter(data, x = "gene1", y = "gene2",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "gene1", ylab = "gene2")

Is this fine or I need to use any log for the scatter plot?

Since Spearman is a rank-based test, it relies on you being able to accurately decide on the ranking of your observations by some metric (usually the magnitude of the numbers). If two observations have identical values (are tied), then they cannot be definitively ranked. Since the ranks are not unique, exact p-values cannot be determined.
e.g. in your data, for gene 2:
                  gene2 rank
TCGA-2Y-A9GS-01A    7.4    1
TCGA-2Y-A9GT-01A    1.0    5
TCGA-2Y-A9GU-01A    2.0   3=
TCGA-2Y-A9GV-01A    0.5    6
TCGA-2Y-A9GW-01A    2.1    2
TCGA-2Y-A9GX-01A    2.0   3=

By way of example, it's possible to do the test with either TCGA-2Y-A9GU-01A or TCGA-2Y-A9GX-01A ranked number 3 - meaning a definitive p-value cannot be determined.
If there are only a few ties, you probably don't need to worry too much. If your FPKMs (or better, TPMs) were calculated with greater precision, then ties would be less likely. Is there a reason you only have these numbers to 1 decimal place?
Log or linear scale doesn't make a huge difference (and neither is 'wrong'). If you choose to plot log counts, you need to consider what happens if you have zeroes in your data.

The presence of a large number of zeros is often the main drawback with using spearmans on gene expression data I find.


text:
Low percentage of reads with consistent barcodes in Split-Seq run
I am following this answer to select reads with allowed barcodes from a Split-Seq run. In particular, I have a whitelist for each round of barcodes (3 rounds) and I filter out reads which have more than 1 Hamming distance from allowed barcodes.
However, after filtering, the number of reads remaining are only the 0.7% of the original number of reads before filtering (169k from 22M). Is this normal, or what may be the underlying issue here?

I would say that it is normal but I'm haven't done splitseq. Could you add some information to help people evaluate it? Like counts distribution across samples, or other plots relevant for QC?

@Llopis I discovered that I was matching one barcode against the wrong whitelist of allowed barcodes. I think this solves the issue, do I need to delete my question?

I don't think you need to delete your question. You could post as an answer how did you find what you where doing wrong, and why it is wrong. If you don't then it might be better to close the question or delete it. But I think it will be useful to know which pitfalls to avoid

In my case the problem was due a wrong set of allowed round1 barcodes, therefore the filtering pipeline was removing most part of the reads. In particular, with the right barcode1 whitelist the total number of reads after filtering is 17M (about 100 times more).


text:
What does it mean for an amino acid to match with itself 99% of the time in a PAM1 matrix?
In this video, the lecturer says:

You see that phenylalanine is probably still phenylalanine because only one percent of amino acids have changed so roughly 99 percent of the time, phenylalanine is still phenylalanine. Now why is it not exactly 99 percent? Because it’s one percent for all 20 amino acids so some are greater than 99 percent and some are less than 99 percent. It averages out to 99 percent.

What on earth does he mean by this? I either do not have the background knowledge or am not thinking hard enough, or I need this to be explained in a way that is more understandable.

It means that the 99% of conservation of amino acids is just a mean value for all the amino acids. That there can be amino acids that are less conserved and some that are more conserved.
So from the 21 basic amino acids, arginine is more conserved than lysine in those alignments. But in general they are conserved the 99% of times.


text:
Seurat Merged objects tSNE - How to paint on original IDs?
I am working with single-cell RNA-seq data, using the R package "Seurat" to cluster and visual data-points.
I had two single cell datasets from which I generated two Seurat objects. I then combined the two using MergeSeurat. I did differential gene expression analysis, performed clustering, and ran a tSNE plot.
Now what I want to do is create that same tSNE plot, but color each cell according to the original object from which it came. How do I do this?

When asking for help, you should include a simple reproducible example with sample input and desired output that can be used to test and verify possible solutions

@MrFlick Don’t you agree that it’s a fairly clear question for domain experts (which, according to your profile, you are)? I’ve flagged it to be moved to bioinformatics.SE.

@KonradRudolph It would be even more clear with an example that I copy copy/paste to test with. Often I help people with packages I've never used before because they include an example that I can play with to get it working. And I don't like answering questions only to have an OP say "this doesn't work with my data"; which happens a lot. So it's just easier to help someone and verify a solution works if a reproducible example is provided. In my opinion there is no question that wouldn't benefit from a reproducible example.

Below a few lines of code that accompany BC Wang's answer.
After using MergeSeurat the sample name will be added to meta data under orig.ident. this can then be used to color the tSNE either using group.by or pt.shape. The former will show colors for each sample, the latter will color each cluster and give sample id another shape.
path1 <- file.path("path to 10x sample 1")
path2 <- file.path("path to 10x sample 2")

S1 <- Read10X(data.dir = path1)
S2 <- Read10X(data.dir = path2)
S1 <- CreateSeuratObject(raw.data = S1, min.cells = 3, project = "S1")
S2 <- CreateSeuratObject(raw.data = S2, min.cells = 3, project = "S2")

s.combined <- MergeSeurat(object1 = S1, object2 = S2, add.cell.id1 = "S1",
            add.cell.id2 = "S2", project = "2samples")

All data analysis goes here:
RunTSNE etc...
Plot the data
TSNEPlot(s.combined, group.by ="orig.ident")
TSNEPlot(s.combined, pt.shape="orig.ident")

If I understand your question correctly, you are trying to present the tSNE plot coloring and labeling the cells representations based on the original sample ID, i.e. which Seurat object it was from.
I think they offered a vignette on that. Please check this link. https://satijalab.org/seurat/merge_vignette.html
Basically, you can add the sample name for the cells in the argument option in MergeSeurat function. Please check the documentation on MergeSeurat function as well.
When I did this, I have a lot of samples to be merged, so I simply saved the sample info for every cell in each sample to different objects. You can also use AddMetadata to do it as well.
Hope that helps.

It is impossible to create the same tSNE plot without knowing which seed you used. tSNE plots are not like PCAs you don't get the same results when using the same data.
But if you want another tSNE plot with a color differentiating from which object each sample comes you'll need a column with these information and color the plot by it.

Seurat allows you to visualize tSNE plots with gene expression data painted on top using the FeaturePlot function. I assume there must be a way to do the same thing with meta data such as original IDs. I just can't find anything online that shows how to do this.


text:
Simulating DNA sequence evolution in R
I hope someone can lend their thoughts on the below code to generate DNA sequences under the Kimura-2-Parameter model of DNA substitution.
The issue is that each time the code is run and the haplotype distribution is examined, there always is a very skewed distribution.
What I would like is for a different distribution to be generated each time the code is run. that is, there should be multiple specimens sharing the same haplotype (instead of most sharing the most common haplotype).
Here is my code:
library(pegas)
set.seed(17)

sim.seqs <- TRUE
length.seqs <- 500
num.seqs <- 100 # number of DNA sequences
subst.model <- "K80" # nucleotide substitution model
transi.rate <- 1e-4 # transition rate
transv.rate <- transi.rate / 2 # transversion rate

if (sim.seqs == TRUE) {

  nucl <- as.DNAbin(c('a','c','g','t'))
  res <- sample(nucl, size = length.seqs, replace = TRUE, prob = rep(0.25, 4))

 if (subst.model == "K80") {
    transi.set <- list('a' = as.DNAbin('g'),
                       'c' = as.DNAbin('t'),
                       'g' = as.DNAbin('a'),
                       't' = as.DNAbin('c'))

    transv.set <- list('a' = as.DNAbin(c('c', 't')),
                       'c' = as.DNAbin(c('a', 'g')),
                       'g' = as.DNAbin(c('c', 't')),
                       't' = as.DNAbin(c('a', 'g')))

    transi <- function(res) {
      unlist(transi.set[as.character(res)])
    }

    transv <- function(res) {
      sapply(transv.set[as.character(res)], sample, 1)
    }

    duplicate.seq <- function(res) {
      num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate) # total number of transitions

      if (num.transi > 0) {
        idx <- sample(length.seqs, size = num.transi, replace = FALSE)
        res[idx] <- transi(res[idx])
      }

      num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate) # total number of transversions

      if (num.transv > 0) {
        idx <- sample(length.seqs, size = num.transv, replace = FALSE)
        res[idx] <- transv(res[idx])
      }

      res
      }

    }

  res <- matrix(replicate(num.seqs, duplicate.seq(res)), byrow = TRUE, nrow = num.seqs)

  class(res) <- "DNAbin"

  # write.dna(res, file = "seqs.fas", format = "fasta")

  h <- sort(haplotype(res), decreasing = TRUE, what = "frequencies")
  rownames(h) <- 1:nrow(h)

}

Output
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

 1  2  3  4  5
96  1  1  1  1

If the code is run multiple times (without the seed), the same pattern emerges: h is always skewed toward the most dominant haplotype.
For example, I want to be able to run the code (without seed) and get something like
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

  1    2    3    4   5
  35  25   20   15   5

i think I would have to specify a distribution for the haplotypes, but the number of haplotypes generated by the mutation process is not known beforehand.
Any thoughts?

Similar to what gringer said: Perhaps I understood wrong your model, but why do you only expect one transition and transversion? Shouldn't that be rbinom( n = 500, ...? So, in a sequence of 500 you expect to know how many substitutions are in all the sequence, not if there is a mutation on the sequence.

Hi Jarrett, the empty double line spacing unfortunately makes your code harder to read than necessary. It would help if you used empty lines only to logically separate blocks, not between all lines.

@KonradRudolph Nice catch! I have fixed the line spacing.

@Llopis The number of substitutions per sequence will not be known in advance. This is determined with rbinom(...).

What are you expecting these lines to do?
num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate)

Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.
I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.
After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.
I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from any existing haplotype) will generally create something new.

These lines compute the number of transitions and transversions across the sequence alignment of a given length according to a binomial distribution. You're correct that more haplotypes will likely be generated when length.seqs and/or transi.rate are increased. However, I'm looking to ensure that there are more sequences per haplotype (so that the distribution is not so skewed).

You're complaining about a skew towards the primary haplotype, i.e. that there are not enough mutations. My recommended fix is to either increase the sequence length or increase the transition rate.

Been playing with the code... still get the same effect when mutation rate and/or sequence length is increased. There should be a better way - I'll ponder it and experiment some more.

My specific application is in DNA barcoding, which itself relies on short (~500 - ~650 bp) DNA fragments. Increasing length.seqs is not really an option - unless, should it work, I then trim the alignment evenly on both 5' and 3' ends to the required length. I'll ponder this.


text:
Extracting sequences from FASTA beginning with common 5' end
I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.
Example of the FASTA file:
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

And what I would like to get from the code is such:
>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.

What is the common region that you are expecting to see?  Also do you want to allow for insertions and deletions or just base edits?

Here is a 'quick' python3 solution, requiring installation of BLASTN and a pip3 install biopython.  The results are output in FASTA format to stdout.
Save the following python3 code in hit_finder.py
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from Bio.Blast.Applications import NcbiblastnCommandline

cline = NcbiblastnCommandline(query=sys.argv[1],
                              subject=sys.argv[2],
                              outfmt="'6 qseqid length qstart'",
                              max_hsps=1,
                              word_size=4)

result = [j[0] for j in [i.split('\t') for i in
                         cline()[0].rstrip().split('\n')]
          if int(j[1]) >= int(sys.argv[3]) and int(j[2]) <= int(sys.argv[4])]

for record in SeqIO.parse(sys.argv[1], 'fasta'):
    if record.id in result:
        print(record.format('fasta'))

Make it executable: chmod +x hit_finder.py.
This script takes four arguments (which in future versions could be enhanced using argparse instead of sys.argv):
arg1 = query_fasta_file (file to search in)
arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
arg4 = start_of_hit_at_5prime (where should the alignment hit start?)

Let's say you wanted to search for AATCTGTAGTACGACTCACTATAGC in your query seqeuences.  Save this sequence in subject.fa as:
>subject
AATCTGTAGTACGACTCACTATAGC

Your query sequences are stored in queries.fa as:
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:
./hit_finder.py queries.fa subject.fa 23 1

To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:
./hit_finder.py queries.fa subject.fa 30 10

To save your results, output stdout to file, e.g.,:
./hit_finder.py queries.fa subject.fa 30 10 > results.fa

This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out.


text:
IonTorrent variant caller failing, lack of ZM: tags in BAM files
I have a lot of IonTorrent sequencing data in BAM format that I want to call variants on, using the IonTorrent variant caller. However, for some of the BAM files (for about 12 samples out of ~150), the variant caller stops due to what it states is lack of ZM: tags in a few individual reads.
The thing is, most of the reads in the BAM file do not have ZM: tags (when viewed through samtools view), but these do not apparently cause problems.
What settings do I need to change to get these BAM files to process, or do I just filter out the problematic reads and analyse the remainder?

It turns out, after speaking to tech support, that these are old files from a previous version of the software so won't run through variant caller.
The files themselves need reanalysing from the initial FASTQ -> BAM with the newest version.


text:
How can I read multiple different regions from a BAM file in R?
I’m trying to process BAM records (of a coordinate-sorted, indexed BAM file) successively in fixed-size, non-overlapping, sliding genome coordinate windows. Unfortunately I am unable to read records via scanBam in different windows (which parameter): the first call to scanBam succeeds but subsequent calls return zero-length results.
By contrast, if I don’t specify a coordinate window and instead use yieldSize when opening the BAM file, I can successfully read records in a loop:
library(Rsamtools)
bamfile = system.file('extdata', 'ex1.bam', package = 'Rsamtools')
process = function (reads) message(length(reads$seq))

bam = open(BamFile(bamfile, yieldSize = 1000L))
param = ScanBamParam(what = 'seq')

repeat {
    reads = scanBam(bam, param = param)[[1L]]
    if (length(reads$seq) == 0L) break
    process(reads)
}

… where process performs the read processing (unimportant for the example).
Now I’m trying to use the following code to work in genomic coordinate windows instead, and it fails as explained above:
windowsize = 1000L

bam = open(BamFile(bamfile))
chr = seqlevels(bam)[1L] # Only do one chromosome for now
chrlengths = seqlengths(bam)

result = lapply(seq(1L, chrlengths[chr] - windowsize, by = windowsize), function (pos) {
    which = setNames(RangesList(IRanges(pos, pos + windowsize)), chr)
    param = ScanBamParam(what = 'seq', which = which)
    reads = scanBam(bam, param = param)[[1L]]
    process(reads)
})

Conversely, it works if I close and re-open the BAM file in each iteration (i.e. before each scanBam call) — but that’s undesirable: it carries a substantial overhead (running on an empty BAM file with 16 chromosomes in its header takes several hours).
It's only returning data from the first iteration, not any subsequent iterations, unless I close the BAM file between scanBam calls.
Since the process function is nontrivial and written in R, I’m stuck with R for now. And there’s no reason why this shouldn’t work in R.
Am I overlooking something? Is this a bug in {Rsamtools}? Is it “by design”?

@gringer I’ve noticed that you’ve added the tag into the question title — why? On Stack Overflow there’s a broad consensus against redundant mentions of the tag in the question title, and in fact I frequently (several times a day, on average) edit to remove such mentions from titles on Stack Overflow (most useless taggings are more egregious, but removing “… in X” is a fair share of the edits I make).

Following the same logic, why not remove BAM from the title? I don't think that "in R" is redundant here. In this case I considered that there were multiple valid Bioinformatics questions that could begin with "How can I read multiple different regions from a BAM file", even with the "R" tag. In this case it seems like you specifically want an R solution, which is why I added it in.

Are you working under linux? and do you absolutely have to do it in R?

@user324810 I mean, this question is three years old. And at any rate, regardless of whether somebody “absolutely” has to do this in R, it should work in R. But for posterity: yes, I absolutely had to do this in R since the process placeholder in my case was a piece of nontrivial logic developed in R over the span of months (in fact, the final implementation was in C++ anyway, the R code was my research prototype for faster iteration). I ended up processing fixed-sized chunks rather than scanning over genomic coordinate windows.


text:
Seurat for clustering bulk RNA-seq?
Is it ever ok to use Seurat for clustering bulk samples?
I am looking at FPKM data from ~750 bulk RNA-seq samples generated using Cufflinks. As suggested for FPKM data, I manually input log transformed data to the @data slot [cd138_bm@data <- log(cd138_bm@raw.data + 1)] and skip the NormalizeData() function. I then use functions FindVariableGenes, ScaleData, RunPCA, FindClusters, RunTSNE, FindAllMarkers in their usual ways to find clusters & cluster markers. My clustering results are quite reasonable, and reflect published work clustering similar samples.
What are the potential pitfalls of using these Seurat functions on bulk data? In FindAllMarkers, would you recommend I use the "negbinom" test? (currently using wilcox) Any other arguments you would recommend changing from the default?

Good question. I removed scrnaseq tag since you have bulk rna sequencing, if I misunderstood somehting, feel free to [edit] it back.

I also want to try this bulk analysis with seurat and I am interested in your methods.
Did you publish this analysis already ?
If it's fine, I would like to ask you get the papers you wrote !

I didn't publish, but I did use seurat on bulk and it worked great.

It's a good question. I processed the sc-RNAseq data and generated some new markers based on gene expression value. In the subsequent analysis, I want to find differential markers in multiple clusters of my dataset based on seurat package.

I'm not sure Seurat is the best tool for this as it was developed for single cell RNA seq data and there are a few intricacies of that type of data that are very different from bulk RNA seq.
For bulk there are really good packages available and corresponding workflows, e.g. limma, edgeR and DESeq2.
The main problems with Seurat for bulk RNA-seq:

Seurat expects counts as input - FPKM are not counts nor are they log counts or log norm counts. It's a specific normalisation method that takes into account gene length and library size and breaks the link between gene counts and variance.
Because of the distribution of reads from a single cell that you get from the sequencer Seurat scales all unique reads (i.e. UMIs) to 10,000 per cell and usually regresses out the nUMI count. Again, becuase you've got FPKMs this makes no sense.
The other functions you can use on bulk data, i.e. tSNE, finding clusters and markers (which is essentially DGE analysis). The one thing to note for the latter is that again because you've got FPKMs they don't follow the negative binomial distribution. Your best bet is limma trend but I'm not sure it's included in Seurat.

You don't need to use Seurat's normalization. You can provide normalized counts or FPKMs and skip the normalization step. See this earlier discussion: https://bioinformatics.stackexchange.com/questions/5115/seurat-with-normalized-count-matrix

That's true but if you look at the original post, OP uses ScaleData, which scales already scaled FPKM values. Even though feeding bulk RNA seq FPKM values to a Seurat@data slot is essentially doable, you are still relying on the downstream processing steps that were designed for log-normalised count data to be meaningful for FPKM values and it's certainly not the case for the FindAllMarkers function.

There are a few issues on the Seurat GitHub that discuss importing TPM/FPKM values. The suggestions are always to just skip normalization and proceed normally. I haven't seen any concerns raised about FindAllMarkers. I don't mean to imply those concerns are not legitimate, but it's interesting they are not raised.

FYI - heard from the author of Seurat, see this github issue
His suggestion to use DESeq won't work in my case, because DESeq requires counts data which I don't have. But it is the green light to generally use the package on bulk data. I am using test.use = "wilcox" in FindMarkers() because it's non-parametric...

I would actually suggest limma-trend (not voom). Compared to a wilcox you'll gain a bit of statistical power due to the empirical bayes moderated t-test.


text:
Transcript Coordinate Ranges to Genomic Coordinates
I have 2 GFF3 files:

Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
Features using chromosome IDs as the landmarks. i.e. "exon" feature types using coordinates from chromosome space.

I would like to transform the coordinates of the features in file 1 into the coordinate space of file 2. i.e. transcript-based coordinates into genomic-based coodinates.
Here is an example for file1:
cat transcript_orfs.gff3
##gff-version 3
##sequence-region Tx.1 1 4000
Tx.1    ORF_finder  gene    1   4000    .   +   .   ID=1
Tx.1    ORF_finder  mRNA    1   4000    .   +   .   ID=2;Parent=1
Tx.1    ORF_finder  exon    1501    2500    .   +   0   ID=3;Parent=2
Tx.1    ORF_finder  CDS 1501    2500    .   +   0   ID=4;Parent=2

gt sketch -addintrons transcript_orfs.png transcript_orfs.gff3

Here is an example for file1:
cat genome.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 genome_assembly gene    1001    6000    .   +   .   ID=Tx
chr3A_part1 genome_assembly mRNA    1001    6000    .   +   .   ID=Tx.1;Parent=Tx
chr3A_part1 genome_assembly exon    1001    3000    .   +   .   ID=Tx.1.exon1;Parent=Tx.1
chr3A_part1 genome_assembly five_prime_UTR  1001    2000    .   +   .   ID=Tx.1.utr5;Parent=Tx.1
chr3A_part1 genome_assembly CDS 2001    3000    .   +   0   ID=Tx.1.cds1;Parent=Tx.1
chr3A_part1 genome_assembly exon    4001    6000    .   +   .   ID=Tx.1.exon2;Parent=Tx.1
chr3A_part1 genome_assembly CDS 4001    5000    .   +   2   ID=Tx.1.cds2;Parent=Tx.1
chr3A_part1 genome_assembly three_prime_UTR 5001    6000    .   +   .   ID=Tx.1.utr3;Parent=Tx.1

gt sketch -addintrons genome.png genome.gff3

I want to convert the coordinates of the features from file1 into chromosome genomic coordinates. I expect to get something like the following output:
cat output.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 ORF_finder  gene    1001    6000    .   +   .   ID=1
chr3A_part1 ORF_finder  mRNA    1001    6000    .   +   .   ID=2;Parent=1
chr3A_part1 ORF_finder  exon    2501    3000    .   +   0   ID=3.1;Parent=2
chr3A_part1 ORF_finder  CDS 2501    3000    .   +   0   ID=4.1;Parent=2
chr3A_part1 ORF_finder  exon    4001    4500    .   +   0   ID=3.2;Parent=2
chr3A_part1 ORF_finder  CDS 4001    4500    .   +   1   ID=4.2;Parent=2

gt sketch -addintrons output.png output.gff3

I have looked into using mapFromTranscripts() from Bioconductor's GenomicRanges library but I have made little progress trying to decipher the manual.

This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.
R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)
#!/usr/bin/env python3

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')

The linked gist contains a class I wrote for doing this on GTF files:
https://gist.github.com/IanSudery/d8349c22823a475ceb489c3e8aeb448e
It uses the GTF class from cgat, which can be found here:
https://github.com/cgat-developers/cgat-apps
You would use it to do this task like so:
from cgat import gtf
from cgat import iotools
import TranscriptCoordsInterconverter

trans_coord_file = gtf.iterator(iotools.open_file(args[1])
transcripts = dict()

for transcript in gtf.transcript_iterator(trans_coords_file):
    transcript_id = transcript[0].transcript_id
    transcripts[transcript_id] = transcript

trans_coord_file.close()

genome_coord_file = gtf.iterator(iotools.open_file(args[2]))
for transcript in gtf.transcript_iterator(genome_coord_file):

    transcript_id = transcript[0].transcript_id
    contig = transcript[0].contig
    strand = transcript[0].strand

    if transcript_id not in transcripts:
        continue

    converter = TranscriptCoordsInterconverter(transcript)
    for exon in transcripts[transcript_id]:
        genome_intervals = converter.transcript_interval2genome_intervals((exon.start,exon.end))
        for interval in genome_intervals:
            new_entry = GTF.Entry().fromGTF(exon)
            new_entry.contig = contig
            new_entry.strand = strand
            new_entry.start = interval[0]
            new_entry.end = interval[1]
            print (str(new_entry)+"\n")


text:
chromPlot R package: remove scales
I am using chromPlot R package to plot some SNPs and their associated score, however, I cannot find a way to remove the scale of such score (circled in the attached picture). Does anybody know how to do it?
Thank you so much.

Please [edit] your question and include the code you use to generate the image.

Setting the legChrom argument to NA will omit plotting a legend.
Found in the manual for chromPlot on Bioconductor for release 3.7.


text:
Compare clusters from different datasets
I have two scRNAseq datasets, and I select a cluster from one of them, and another cluster from the other. I want to find differentially expressed genes between the two clusters. What would be the right way of doing that? Right now I am using DESeq2 in the following way. I select only the selected clusters' cells, merge them together forming a matrix and then feed them into DESeq2:
mtx <- merge(data_1,data_2,by="row.names",all=TRUE)
mtx[is.na(mtx)] <- 0
mtx_conv <- matrix(as.numeric(unlist(mtx[,2:ncol(mtx)])),nrow=nrow(mtx))
mtx_int <- round(mtx_conv, 0)
rownames(mtx_int) <- as.character(mtx$Row.names)

cluster <- as.factor(c(rep(0, dim(data_1)[2]), rep(1, dim(data_2)[2])))
colData <- data.frame(
  Cluster = cluster
)

dds <- DESeqDataSetFromMatrix(countData = mtx_int,
                          colData = colData,
                          design = ~ Cluster)
dds <- DESeq(dds)

I am not sure though that it is the correct way to go because I am not using the full original datasets. The reason is that I want the computations to run faster. Another way of doing that is probably combining the whole two datasets, defining 3 clusters: two groups I am interested in, and all other cells, and then running the analysis. Would it differ significantly? Am I on the right track at all?
I have another idea: to use Seurat package. There, I would run CCA algorithm to align two full datasets, and then run FindMarkers function between the two clusters. I feel like it may be wrong, because the two datasets may need to be re-normalized together but Seurat does not seem to be doing that:
https://satijalab.org/seurat/Seurat_AlignmentTutorial.html
Any suggestions would be greatly appreciated.

I've been working with the multi-CCA alignment workflow in the Seurat package and am very impressed with functionality. As usual, the CCA tutorial is excellent for an alignment workflow and will get you off the ground quickly. You can even use the DESeq2 algorithm for determining differential expression from within Seurat.
When you use the function RunCCA() you can pass rescale.groups=TRUE if you wish to re-scale.
Abstracted workflow

Data preprocessing and gene selection
Define a shared correlation space with canonical correlation analysis
Remove dataset specific cells
Align correlated subspaces using dynamic time warping
Integrated analysis across datasets (clustering, trajectory building, etc ...)

Butler, A., and Satija, R. (2017). Integrated analysis of single cell transcriptomic data across conditions, technologies, and species.