text:
RIP-seq analysis?
Given an experiment consisting of an input (baseline RNA) and IP (pulldown to find RNAs bound to certain protein of interest)... Is a DE analysis performed over the RNA-seq data from the samples (lets say with EdgeR or DESEQ2) suitable to reveal the preferentially bound RNAs? What other software tools would you recommend?

You need to be careful of terminology. To me, a RIP-seq experiment involves a pull down, followed by a RNAseq library prep. Thus, the whole transcript is captured, not just the binding site of the protein (as in CLIP-Seq, HITS-CLIP, PAR-CLIP, iCLIP or eCLIP). Thus "peak-callers" whether they be designed for calling protein-DNA or protein-RNA binding sites are not suitable as the signal will not be punctuate or peaky. Of the methods mentioned by @Devon Ryan, only RIPSeeker seems setup to deal with this sort of data.
In some of the papers that talk about calling RNA-protein interactions they lump RIP-seq in with the CLIP techniques (see piranha and ASPeak papers referenced by @Devon Ryan). They appear to be talking about a technique where the RNA is fragmented before it is pulled down. Thus you would capture peaky binding sites (actually, the protocol would be remarkably similar to CLIP).
You should work out which of these you have.
As for using DESeq/EdgeR etc, in the RIPSeeker paper they say:

Furthermore, programs for de novo transcript assembly followed by
  differential expression (DE) analysis, such as the Cufflinks/Cuffdiff
  suite (15,16), and for DE on a set of known transcripts, such as DESeq
  (17), may appear applicable to RIP-seq analysis. Unlike peak-calling
  strategy, however, the transcript-based methods assume the full
  transcriptome being sequenced at a fairly deep coverage (as usually
  the case in RNA-seq) and thus may be sensitive to background noise
  typical to the IP-based protocols, which is due to both the
  non-specific RNA interactions with a protein of interest and the
  non-specific RNA input from the pull-down of the (mutant) control
  (Supplementary Figures S3 and S4).

I don't know to what extent RIPSeeker performs better than the naive DESeq approach as its not included as one of the comparators in the paper.

How important is peakiness though? Macs2 seems to deal well with broad peaks. Couldn’t it be used here?

There are number of reasons not to use a protein-DNA peak caller, many detailed in the RIPSeeker paper. Many of these are to do with splicing: fragment sizes estimated from mapped RIPseq will not reflect genuine fragment sizes, using either extended reads or fragment midpoints will likely leave you in the middle of an intron.
MACS2 broadpeak works by merging narrow peaks connected by sub-significant sequence. Its not clear this is a good description of the data you'd get in RIP. The RIPSkeeker paper compares to Macs, and looks like it out-performs by quite a way to me.

You can try doing standard differential expression, but I worry that the between-sample normalization will work poorly. Personally, I would do peak calling instead, followed by diffBind. You have a few tools to choose from when it comes to this. In the past, I've rolled my own methods for this using MACS2 and genomic alignments (I then converted those to bedGraph files where entries are transcripts, so peak calling is then in transcript-coordinates).
It's probably easier, though, to use premade software. There are a number of packages out there for RIP-seq. These include piranha, RIPSeeker and ASPeak. Particularly if you're just starting out, you'd be well served with using one of those tools.

I've got the RIP-seq analysis outputs. We used an isotype control antibody VS antibody against an RBP to pull down corresponding RNAs!
I did the following pipeline: Hisat2 (mapping the reads after QC) -> htseq-count (getting read counts) -> DESeq2 (getting enriched RNAs in test VS isotype control) -> volcano plot shows up-regulated and down-regulated genes in the test sample.
Logically, I should have seen just enriched RNAs (up-regulated genes), but surprisingly there are some RNAs down-regulated in the test sample (or up-regulated in the isotype control sample)!!!!!
How can we explain and interpret this???? it seems some genes are pull down in isotype control sample as well. Are these background? Is my computational pipeline right?
BTW, it's worth to mention that all enriched RNA (called up-regulated genes in DESeq2 pipeline) are significantly relevant and meaningful to my experiment. However, I am not sure how much confident I can be on this output? And how should I explain down-regulated RNAs? because normally we just expect enriched RNAs in RIP-seq experiment.
Any thoughts?


text:
Extract LINE-1 methylation status from WGBS data
I've been trying to figure out the correct way of extracting methylation status of particular elements from WGBS data, and found several possible approaches, but would like to discuss which one could be the best:

using software dedicated to the extraction of these elements on WGS (like TraFiC or MELT and then quantify methylation (like Bismark)
using bedtools to scan bam file with a filtered version of UCSC RepeatMasker table in bed format) and then quantify methylation
I've seen the use of custom perl scripts on fastq files like in this paper, in which the Supplemental Methods Section states: "A custom Perl script was written to assay methylation of CpGs within a previously published motif of LINE-1 elements (Yang et al., 2004). Methylation levels for each of the four CpGs within the consensus LINE-1 element were determined globally for each WGBS sample prior to alignment using Line1_FASTQ.pl. Dup15q and control were compared and tested for significance using an unpaired t-test."

Any thoughts on this?

Have you checked out WGBS Tools on github? All of the custom scripts used in the paper that you referenced are in that repository. Each tool has a brief description of how to use it.

Via BEDOPS:
$ bedmap —-echo —-echo-map-id-uniq <(bam2bed < reads.bam) <(rmsk2bed < repeats.rmsk) > answer.bed

After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:
element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.

source: Stack Exchange

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

sort: votes


text:
How to convert data in gmt format to dataframe?
I downloaded c5: gene ontology gene sets file from http://software.broadinstitute.org/gsea/downloads.jsp
I opened the "c5.all.v6.2.symbols.gmt" file in csv format and It looks like below:

I want to convert the .gmt file into a dataframe. And it should look like below:

can anyone say how I can do that. Thanq

This might be an XY problem (you are asking about a problem, but you really have another one). Why do you need this data in data.frame format? Also consider that the GO data is a directed acyclic graph, which has a special structure that should be taken into account

for some analysis I want the data to be looked like the output I posted.

Also note that the GO annotations imported by MSigDB are quite old. If you are interested in associations between GO and genes (updated) you should use other sources

other resources like?

Like the GO database directly or PANTHER. Or to retrieve the results biomart or its package biomaRt

Ok. Answering my own question.
I did this way. And I got the output I want.
install.packages("msigdbr")
library(msigdbr)

m_df = msigdbr(species = "Homo sapiens", category = "C5")
head(m_df)

Output:
gs_name gs_id gs_cat gs_subcat human_gene_symb~ species_name entrez_gene
  <chr>   <chr> <chr>  <chr>     <chr>            <chr>              <int>
1 GO_14_~ M184~ C5     MF        AANAT            Homo sapiens          15
2 GO_14_~ M184~ C5     MF        AKT1             Homo sapiens         207
3 GO_14_~ M184~ C5     MF        ARRB2            Homo sapiens         409
4 GO_14_~ M184~ C5     MF        BAD              Homo sapiens         572
5 GO_14_~ M184~ C5     MF        DAB2IP           Homo sapiens      153090
6 GO_14_~ M184~ C5     MF        DDIT4            Homo sapiens       54541
# ... with 2 more variables: gene_symbol <chr>, sources <chr>

found this solution from here https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html

perfect! thanks.

You can import the csv file as a table. This does not require further libraries:
#you're reading a csv file, using tab as field separator and considering the first line as the headers of the data table
data <- read.csv(filename, sep="\t", header=TRUE)

then what you need are just two columns of your data:
colnames(data) #will give you the names of the columns that can be accessed
#to access two columns and merge them into a selectedColumns table
selectedColumns <- cbind(data$colName1, data$colName2)

The question is tagged with R so it is an R data.frame.

you're right! I missed it

You can use GSA.read.gmt function from GSA package. The following code can be used to convert the file to a dataframe. Just ignore the warnings.
Original_response
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
gene_names <- unlist(data$genesets, use.names=FALSE)
your_dataframe <- cbind(data$geneset.names,gene_names)
colnames(your_dataframe) <- c("Pathways","Genes")
mydataframe <- as.data.frame(your_dataframe)
head(mydataframe)
                                               Pathway   Gene
1            GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2                           GO_CARDIAC_CHAMBER_DEVELOPMENT POLR2J
3 GO_DNA_DEPENDENT_DNA_REPLICATION_MAINTENANCE_OF_FIDELITY  CTDP1
4                                      GO_CIRCADIAN_RHYTHM   RDBP
5              GO_PHOSPHATIDYLSERINE_ACYL_CHAIN_REMODELING COBRA1
6                               GO_SPINAL_CORD_DEVELOPMENT   RSF1

Edited_response
You can use the code given below to achieve your desired output:
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
len_vec=c()           # Now create a vector for containing the length of genes at each position
len_vec[1] = 3
for(i in 1:length(data$genesets)){len_vec[i] <- c(length(data$genesets[[i]]))}
pathway_vec <- unlist(Vectorize(rep.int)(data$geneset.names, len_vec),use.names = FALSE) # Now create a vector for all the pathways in the data
desired_df <- as.data.frame(cbind(pathway_vec,unlist(data$genesets,use.names = FALSE))) # This gives your desired dataframe
head(desired_df)
                                    Pathway      Gene
1 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2J
3 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION  CTDP1
4 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RDBP
5 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION COBRA1
6 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RSF1

sorry, this is wrong. Please check the output I want in the above question. In your output Genes POLR2J, CTDP1, RDBP, COBRA1, RSF1 shows different pathways.

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


text:
How to convert featureCounts to FPKM?
I have seen many posts regarding counts to RPKM and TPM. I haven't seen any post for counts to FPKM.
I have RNA-Seq data which is paired-end reads. Extracted the counts using featureCounts for all the samples.
There is a function to convert counts to RPKM: using the gene_length
rpkm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(counts) * 1e6
}

I know that RPKM is mainly used for single-end reads data. Do you think I can use the above function for converting counts to FPKM as my data is paired end? [TCGA data: HTSeq counts are converted to HTSeq FPKM]
Somewhere I have seen that the same function applies for single end data which will be RPKM and also for paired end data which will be FPKM. Am Ir right or wrong?
If not can anyone show some function or code to convert counts to FPKM please.
thanq

Why do you need RPKM/FPKM? Haven't you heard that these values are obsolete?

I guess OP wants to follow the most up to date procedures, but seeing that in TCGA they used FPKM is trying to work like them. I thought that FPKM were ok: see this https://bioinformatics.stackexchange.com/a/4270/48, or this https://bioinformatics.stackexchange.com/a/69/48. @beginner I hope the links are helpful, update the question otherwise

@Llopis thanks for the links. In this https://bioinformatics.stackexchange.com/questions/66/how-to-compute-rpkm-in-r/69#69 I see that FPKM function using effective length...it is gene_length right?

@b.nota Yes, but as TCGA data is FPKM I wanted to use FPKM data.

@beginner See this question and answers: https://bioinformatics.stackexchange.com/q/367/48

I have seen many posts regarding counts to RPKM and TPM.

There’s your answer then: FPKM = RPKM. It’s simply a more accurate name.
Speaking of RPKM for paired-end data is discouraged because the reference to “read” in this context lends itself to ambiguity. But mathematically the quantity is the same: we are counting fragments, not individual reads (of which each fragment has two, for paired-end data).
But as mentioned in the comments, there are good reasons against using FPKMs.

thanks for the answer. I have a question. I have Ensembl gene ids as rows and samples as columns with counts. I have an other column gene_length. I'm trying to convert those counts to TPM. In this post [https://gist.github.com/slowkow/6e34ccb4d1311b8fe62e] they used gene_length to convert counts to TPM. But I have seen one of your post in Rpubs [http://rpubs.com/klmr/rnaseq-norm] where you used effective length. Which one is the right one for converting counts to TPM?

In every case, the correct number is the effective length, i.e. the transcript/gene length minus the average read length plus 1 (this is the “effective” length because it represents the number of possible alignment positions of a read to the transcript). That said, in practice it should make very little difference.

I see that in the post it is given EffectiveLength = Length - mean_fragment_length + 1; And how the mean_fragment_length is 300? One more question can I use this for Ensembl gene ids? or only for transcripts?

@beginner The fragment length is determined by your sequencing protocol. It’s (more or less) the distance between the two paired reads (including the read lengths, minus any overlap). For a typical long RNA-seq library prep it’s often in the order of 300–400 bp. You can run Picard’s CollectInsertSizeMetrics tool on your aligned BAM files to find the mean fragment length.

Sure. I will. I'm working with TCGA data. could you please tell about my other question I the previous comment.

You can use this on any annotated feature you like, both genes and transcripts. There are some fundamental issues with using genes (essentially, you lose information about alternative transcripts, and gene-level statistics by necessity need to average across transcripts) but that isn’t specific to quantification, and in practice it’s usually not an issue.

You can use countToFPKM package.
This package provides an easy to use function to convert the read count matrix into FPKM matrix; following the equation in

The fpkm() function requires three inputs to return FPKM as numeric matrix normalized by library size and feature length:

counts A numeric matrix of raw feature counts.
featureLength A    numeric vector with feature lengths that can be
obtained using
biomaRt.
meanFragmentLength A numeric vector with mean fragment lengths,
which can be calculate with Picard
using CollectInsertSizeMetrics.

See https://github.com/AAlhendi1707/countToFPKM

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


text:
Calculating Z-score from logCPM values using edgeR
I have the raw counts for RNA-Seq data. I converted counts data to logCPM using edgeR package.
Lets say I have a dataframe A with 15000 genes as rows and 100 samples as columns with counts data. I have used following code from edgeR tutorial.
logCPM <- cpm(A, prior.count=2, log=TRUE)
head(logCPM)[1:5,1:5]
               Sample1          Sample2          Sample3         Sample4
A1BG          -1.019429        0.8333739        -1.192960       -0.6717287
A2M            7.622280        7.8057767         7.932570        8.6409077
A2ML1          1.004541        1.0909119         5.462613        6.4149348
A4GALT         2.126502        1.8403954         5.322601        2.9469345
AAAS           5.785888        5.1639501         5.259789        5.3248026
               Sample5
A1BG          0.3121975
A2M           9.1949911
A2ML1         9.4740594
A4GALT        2.9214979
AAAS          4.3708063

After converting counts to logCPM I see there is another command for scaling which is used for heatmaps.
logCPM <- t(scale(t(logCPM)))

It is given there that scaling each row (each gene) to have mean zero and standard deviation one in page 15 edgeR
Is this nothing but z-score? or Is there any otherway to calculate z-score from logCPM data?
thanq

Yes this is nothing but a z-score. The scale function calculates column z-scores, but since you use t (transpose), it is actually performed on rows. There are many other ways to calculate z-scores, but if done correctly they will all give the same results.


text:
Question about Co-expression analysis and finding targets for lncRNAs
I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)
I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.
Few doubts:
1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?
2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?
3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.
Thank you

I think you should ask a separate question for your third question (remember to show what you have tried/searched...). For the question 1: which is your biological question? Are you interested in correlation between lncRNAs and mRNAs regardless if there is a tumor or not? Or are you interested in knowing if the co-expression in normal samples are different from the tumor? For question 2: You should be aware that you can find lots of false positive on those correlations. You should at least restrict yourself to those significantly correlated. But I wouldn't recommend to follow this path...

@llrs I'm interested in co-expression in tumor, may be I will do spearately. For the second, ofcourse I go based on pvalue.

In this case it is better if you only use the tumor samples. If you go based on p-value consider correcting the p-value to reduce the amount of false positives

And for WGCNA, after module detection, there are some important genes Im interested in were kept in the grey module which is unassigned genes. What to do in this situation?

You can then say that these genes are not co-expressed following a clear pattern. Either they don't follow a pattern, or are the ones regulating the expression of the other modules, or you don't have enough power/samples or it is too noisy to find a signal... The important genes might be for other reasons not for co-expression between mRNA and lncRNA (or they still are but you can't detect it)

You have several options to approach this with WGCNA (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.
Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).
I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.
I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.

Thanks a lot for the answer Peter. With WGCNA, I detected the modules and from that modules having mRNAs and lncRNAs how to detect hub genes inside modules?

Hub genes are the ones with high intramodular or eigengene based connectivity. Use the function signedKME to calculate kME values (it's really just a correlation of module eigengenes with individual expression profiles) and choose the top genes.

There is no hard and fast criterion for what constitutes a hub gene; I'd probably look at the first maybe 20 genes for a smaller module (less than say 200 genes) and maybe a few more (say up to 50) for larger modules; if you have a module of more than a thousand gene, looking at the top 100 may be useful. You can also go by kME (i.e., correlation); anything higher than 0.7 is pretty strong, and sometimes you can go down to say 0.6 to still find useful answers.

Thanks peter. After module detection, I used this code. Do you think this is right?         MEList = moduleEigengenes(datExpr, colors = dynamicColors);
MEs = MEList$eigengenes;
datKME=signedKME(datExpr, MEs);
head(datKME);
FilterGenes= abs(datKME$kMEbrown)>.8;

table(FilterGenes);
dimnames(data.frame(datExpr))[[2]][FilterGenes];

Hi Peter, Could you please check the above code whether that is right way to get hub genes?

As far as I can tell (the formatting is messed up) it should correctly give you the top hubs in the brown module.

A cursory search to bedtools documentation will reveal the bedtools closest feature - which might be exactly what you are looking in your third question. You can download the list of protein-coding genes from Ensembl biomart or UCSC table browser in bed format, and also convert the list of your lncRNAs into a bed compatible file. Rest is easy:
bedtools closest -a lncRNAs.bed -b proteinCoding.bed [-iu/-id] -D ref

2) If not WGCNA, if I use cor function in R with method Pearson, on
  what cutoff I should select the target genes?

Good call, this is a nice approach. It is a good question because the central issue is that a purists approach of 0.7 results in insufficient information to understand the behaviour of the data.
The guys here were using 0.4, in my experience co-variance of 0.5 and above is good evidence, whilst 0.7+ is strong evidence.
I've done this (not on RNAseq) and I used 0.25 to 0.45 and 0.45 and above then used a scatterplot with one half of the graph showing low co-variance and the other half strong co-variance. So in one sense I ignored my own advice (0.5+). In my defence its easy to reassign the boundaries. Each scatterplot (1000s) was then subject to tSNE analysis to assess the group behaviour. Its a different context but gives an idea.
Networks are cool for visualisation but they are not particularly  quantitive for most (not all) network theory.

Just to explain to everyone else, what happens if you get boundary effects with the data and you are assigning a discrete boundary on a continous distribution (its a negative binomial distribution in my case). I observed a vast amount of results  occurring between 0.2 to 0.25 - but 0.2 isn't good evidence, so that data has to be lost. These sort of boundaries were in days before genomic throughput. Thus in my context I didn't consider a hard boundary particularly informative.


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:
How to downsample some of the samples in RNA-seq data?
I have 40 samples and these are into two groups. I would like to perform a differential analysis between two groups. The library size of the samples is very low. But there are two samples (GroupA_12 and GroupB_2) in which the number of counts is too large compared to other samples. Please check the library sizes below:

I wanted to downsample these two cases and then would like to perform differential analysis. May I know how to downsample those two samples.
Any help is appreciated. thanq.

You have a few options:

Downsample the fastq files and rerun the entire analysis. You can do this with seqtk sample.
Downsample the BAM files, which you can do with samtools view -s.
Divide all of the counts in the counts files by some factor and round that to an integer.

I personally prefer option 2, since it's quick and doesn't usually have any side-effects unless you use two-pass alignment or did the quantification with something like salmon. Option 3 is always tempting, but I have to say I've never been convinced that this consistently produces similar enough results to options 1 or 2.

thanks for the reply @Devon Ryan, Yes, I feel that the 3rd option is much easier. To use the 2nd option, I have to mention that the raw read counts of the transcripts were extracted from the stringtie step using prepDE.py python script [https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual] Will there be any problem?

I don't use stringTie, so I can't comment with certainty.

ok. And if I have to use the 2nd option, Is this right? samtools view -s number GroupA_12.bam >  GroupA_12_new.bam What I should give the fixed number for downsampling?

If you need 20% then 1234.2. 1234 is just a seed for reproducibility, .2 is the fraction to keep.

Small question. How to decide the fraction of reads to keep? I mean 0.2 or 0.5 or something else? On what basis I should give the fraction?

You’re trying to make the library sizes similar.

yes...so any cutoff is fine?

Just eye ball a fraction from that graph.


text:
How to calculate logCPM across all samples?
Using edgeR for differential analysis between Tumor and Normal gave me differential expressed genes with logFC, logCPM, PValue and FDR.
From the details of glmTreat function I see that logCPM is average log2-counts per million, the average taken over all libraries.
The output table I got from glmTreat:
Geneid     logFC  unshrunk.logFC    logCPM  PValue       FDR
Gene1   4.985298235 5.573931983 1.236660486 1.55E-54    1.87E-50
Gene2   4.613438634 5.52126484  0.920343233 2.59E-53    1.56E-49
Gene3   5.250601296 5.356432653 2.751294034 2.48E-50    9.95E-47
Gene4   4.943049159 5.379741165 1.361393757 2.06E-45    6.19E-42
Gene5   6.100121754 6.117580436 5.774385315 2.95E-43    7.11E-40
Gene6   4.722697891 5.320461275 1.120685402 3.05E-42    6.11E-39
Gene7   5.246129853 5.497012001 1.902992053 2.78E-40    4.78E-37
Gene8   3.878773277 4.956636276 0.776208741 8.66E-39    1.30E-35
Gene9   4.441752496 4.94930499  1.132652682 4.37E-38    5.83E-35

To calculate logCPM given in the above table manually I did like this:
logCPM <- cpm(y, prior.count=2,  log=TRUE)

I got the logCPM values for all the genes. Then I calculated Average across all the samples for each gene. I got like below:
Geneid  Average
Gene1   0.686560246
Gene2   0.617115826
Gene3   1.075975225
Gene4   0.692050878
Gene5   1.277556065
Gene6   0.638358189
Gene7   0.689323163
Gene8   0.60700396
Gene9   0.662115092

Why there is so difference? Where I'm doing wrong?

Did you also try aveLogCPM(y, normalized.lib.sizes=TRUE, prior.count=2, dispersion=NULL, ...)?

glmTreat() literally uses aveLogCPM under the hood in fact.

@DevonRyan Yes you are right. thanks. It is calculated across all samples. I want to calculate aveLogCPM for each condition differently (Tumor and Normal separately) How to do that? Just taking the raw counts of Tumor condition and using aveLogCPM function on it is fine?

If you want the values by group then subset y to contain the samples of interest and then feed that to aveLogCPM().

yes I took the matrix of only Tumor samples and genes with raw counts and applied aveLogCPM()

This post on the Bioconductor mailing list explains the difference: https://support.bioconductor.org/p/100463/

Thanks for your response. Could kindly provide a summary of the link here? Link answers are generally not encouraged.

In short, Aaron Lun explains that the output of aveLogCPM is closer to
log2(rowMeans(cpm(y)) (where the log-transformation is done after the averaging), than to:
rowMeans(cpm(y, log=TRUE))


text:
What could be the reason for the samples not clustering?
I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in Stringtie manual. With the counts data I have used edgeR for Differential analysis.
Before that I have seen MDS plot showing in the following way:

I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.

Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?

Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.

Yes, you are right. If some thing went wrong what all should I check anything like mapped reads percentage or anything else?

Mapped percentages, size factors, and perhaps sparsity (percentage of genes with no counts). Basically you want to ensure that they're all similar. Also, try redoing the counts with featureCounts. Then you can eliminate an oddity in the counting as an issue.

Which is the mapping percentage in the above comments? It gives overall alignment rate.

The alignment rate is the mapping percentage.

Counts with featurecounts should I use bam or sorted.bam?

It doesn't matter which you use.

I initially checked on one bam file. Geneid, Chr, Start, End, Strand, Length, /path/to/A.sorted.bam. These are the columns in the featureCounts output. Is the last column read counts? Do you think this is fine? featureCounts -a annotation.gtf -F GTF -p -s 2 -T 8 -P -o countMatrix_all.txt /pathto/*.sorted.bam

Yes, the last column is the counts (I wish there were an option to simplify that output to just geneid and count). You probably don't need the -P option, not that it's likely hurting anything. BTW, that could be its own post here.

One last question. When I did samtools flagstat on one of the bam file i see this. 202671876 + 0 in total (QC-passed reads + QC-failed reads)
128036784 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
179886195 + 0 mapped (88.76% : N/A)
74635092 + 0 paired in sequencing
37317546 + 0 read1
37317546 + 0 read2
45092910 + 0 properly paired (60.42% : N/A)
48583902 + 0 with itself and mate mapped
3265509 + 0 singletons (4.38% : N/A)
621892 + 0 with mate mapped to a different chr
377778 + 0 with mate mapped to a different chr (mapQ>=5) In this what is mapped and properly paired?

Please post things like this as stand alone questions (make sure to google them first, this is answered elsewhere).

The largest variation differences in total expression will be seen in cell population differences. It's a good thing that your conditions aren't clustering on the first and second principal components when all genes are considered, because it means that cell population is probably not a factor in the variation of expression.
There may be some other principal component that does separate the conditions, but you might need to drill down fairly deep into the data. This will depend on how much the treatment conditions contribute to the total expression variation. We got lucky with our Th2 RNASeq study in that cell population was mostly on the first and second principal components, and treatment condition appeared on the third principal component.
If you want to force condition clustering, identify the genes that are differentially expressed, and only use those genes for the principal component analysis.

I'm not sure how useful doing a PCA on DE genes is. You already know your samples differ by group in those genes, after all.

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


text:
5' and 3' bias in Rna-seq data
I'm working with rna-seq samples. I see 5' bias and also 3' bias in the per-base sequence content plot. From this link I see that the bias at the start of the sequences appears to be the result of biased selection of fragments from the library. And this won't effect any downstream analysis. But In all of my samples I also see bias at the end of the sequences. What can be the problem? Samples are strand-specific sequenced with poly-A selection protocol.
In the Transcript coverage profile I see 5'-3' bias for all my samples is: 0.9
Showing the per-base sequence plot for one of the sample: For all the samples I see that 5' and 3' bias.
When I looked at Adapter content this is what I see:
No samples found with any adapter contamination > 0.1%

I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.
It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was done or there was a bit of degradation. The primary thing to ensure is that the bias is similar across samples/groups. If that's the case then you don't have to worry that your downstream analyses are going be affected by some sort of per-group bias.

I'm using hisat2 and yes I see the bias is similar across samples. Thanks for the answer.

Ah, for hisat2 you might get a bit better results if you trim off the last 1-2 bases. This is because hisat2 does global alignments (i.e., it doesn't do soft-clipping). Or at least that was the case last time I checked.

Cool. I will consider this and will trim few bases in the 3' end. Does it need any trimming in the 5' end also?

Nope, the 5’ end should be ok.

I thought tha hisat2 had a --local option, but when I run hisat2 -h | grep local, I only get --ma <int>         match bonus (0 for --end-to-end, 2 for --local). No --local proper. The help does not look coherent.

@bli The joys of bioinformatics package documentation...

I know this is an old post, but if this plot is from after trimming I would suggest a different explanation: some trimming tools remove poly-A sequences from reads. If that's the case then any read ending with A will have that removed, this leads to a 0% A base content in the final base position (and a dropping %A in the final couple of bases).
This should have no effect on alignments or analysis.


text:
Greedy Motif Search Using Probability Matrices
In my greedy search with pseudocounts algorithm in my bioinformatics course, I did not follow the pseudocode since I wanted to solve the problem in my own way. Unfortunately, although my answer is correct, the answer we are supposed to give is the first correct answer, not the last.
The problem works as such: Given t strings, you make a profile matrix from the first k-mer in the first string, and then you look at the next string and using that profile matrix, you find the best k-mer match, and then you combine that k-mer with the one from the first string to make a new profile matrix, etc. You eventually do this for all k-mer windows in your first string of DNA. The way I did it was by making a probability profile matrix from the best k-mer groups in each case, and then multiplying the greatest probabilities in each column. This gives me a score, and if the next score is higher (thus, more probable) that group of k-mers is used. Below is my code:
 def highest_probability(string, k, matrix):
    score = 1
    temp_string = ""
    c = 0
    best_string = ""
    temp = 0
    while c+k < len(string):
        for i in range(c, c+k):
            if string[i] == 'A':
                score *= matrix[i-c]['A']
                temp_string += 'A'
            elif string[i] == 'C':
                score *= matrix[i-c]['C']
                temp_string += 'C'
            elif string[i] == 'G':
                score *= matrix[i-c]['G']
                temp_string += 'G'
            elif string[i] == 'T':
                score *= matrix[i-c]['T']
                temp_string += 'T'
            if i == c+k-1:
                if score >= temp:
                    temp = score
                    best_string = temp_string
        temp_string=""
        score=1
        c+=1
    print(best_string)
    return best_string

def GreedyMotifSearch(DNA, k, t):
    best_motifs = [i[0:k] for i in DNA]
    score=0
    c=0
    while c+k<len(DNA[0]):
        k_mer = DNA[0][c:c+k]
        motif=[k_mer]
        matrix = make_matrix(k_mer, matrix=[], i=0)
        for i in range(1, t):
            compare_string = DNA[i]
            good_motif = highest_probability(compare_string, k, matrix)
            motif.append(good_motif)
            matrix = make_matrix(good_motif, matrix, i)
        temp_score=1
        for dictionary in matrix:
            temp_score *= max(dictionary.values())
        print(temp_score)
        if temp_score > score:
            score = temp_score
            best_motifs = motif
        c+=1
    return best_motifs

def make_matrix(string, matrix, i):
    if len(matrix)==0:
        for ch in string:
            new_dict = {'A':0, 'C':0, 'G':0, 'T':0}
            new_dict[ch] = 1
            matrix += [new_dict]
        print(matrix)
        return matrix
    elif len(matrix)==3:
        for j in range(len(string)):
            matrix[j][string[j]] += 1
        Sum = [sum(matrix[0].values()),sum(matrix[1].values()),sum(matrix[2].values())]
        for columns in range(len(matrix)):
            for keys in matrix[columns]:
                matrix[columns][keys] /= Sum[columns]
        print(matrix)
        return matrix

DNA=['GGCGTTCAGGCA',
'AAGAATCAGTCA',
'CAAGGAGTTCGC',
'CACGTCAATCAC',
'CAATAATATTCG']
k=3
t=5
print(GreedyMotifSearch(DNA, k, t))

The output gives me ['TTC', 'GTC', 'TTC', 'GTC', 'TTC'] which is equally correct to the answer given in the problem output which is TTC,ATC,TTC,ATC,TTC. Unfortunately, my answer is the last correct occurrence. I tried adjusting my highest_probability function with if score is greater than temp instead of greater or equal, but that gives me the answer ['GGC', '', '', '', ''] which I don't understand at all. My question is 1) Why is this the case? and 2) How can I fix this? And does anyone have a simpler way of doing this?

Assuming you're talking about Implement GreedyMotifSearch from Rosalind, while I haven't actually finished this particular one, I've got an idea why TTC,ATC,TTC,ATC,TTC doesn't work. Here's my explanation from the "Questions" section of that problem:

The profile is created from the first kmers in each string, and the most likely string from those is CAG. The answers will generally be similar to this.

I've been trying to do Rosalind by only using the information provided in the questions; this is one of the many problems where there wasn't enough information in the Question for me to arrive at a correct solution when I tried it. It's possible that the debug datasets were added to this question since I last looked at it, as they seem to describe this situation:
TEST DATASET 1:
Input:
3 4
GCCCAA
GGCCTG
AACCTA
TTCCTT

Output:
GCC
GCC
AAC
TTC

This dataset checks that your code always picks the first-occurring Profile-most Probable k-mer in a given sequence of Dna. In the first sequence (“GCCCAA”), “GCC” and “CCA” are both Profile-most Probable k-mers. However, you must return “GCC” since it occurs earlier than “CCA”. Thus, if the first sequence of your output is “CCA”, this test case fails your code.


text:
How to convert DNAbin to FASTA in R?
I am trying to convert DNAbin files to fasta format. Rationale: I want to use the fasta file to calculate non-synonymous/synonymous mutation rate.
my_dnabin1 is a DNAbin file of 55 samples and I am using the following code to convert it into a fasta file.
dnabin_to_fasta <- lapply(my_dnabin1, function(x) as.character(x[1:length(x)]))
This generates a list of 55 samples which looks like:
$SS.11.01
[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"
$SS.11.02
[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"
and so on...
However, I want a fasta formatted file as the output that may look something like:
>SS.11.01 ttacctga
>SS.11.02 ttacctga

The following should work:
cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

A bit of explanation, starting from the middle outwards:

Collapse the letter arrays inside dnabin_to_fasta into single strings
Append ">" to the front of the names of each sequence
Concatenate the results of steps 1 & 2 together, separating with a line break
Write each concatenated string to a file, separating by a line break

Thank you so much! It solved a major part of my problem. When I print it to the screen, I get the exact output I want. However, it does not produce any output files. Could you please help me with this, too?

Note that R will save files in the "working directory", which is not necessarily the directory that the last file was loaded from. If you're using cat with the file option, it's probably saving somewhere. You can find out where by using the getwd() command.

Thank you so much! It solved my problem! I was stuck in it since a long time.


text:
Read and write FASTA files with more information than id and sequence
I am trying to read a fasta file, manipulate is in Python (using BioPython) and then write it back.
The format of my sequences is like:
>k119_5 flag=0 multi=141.0706 len=473
AGGTTAGTCAGCACCGTTTCCGTGGTGCTGCCTTTCGCTTCAAAACCGACGGCGTCTATTACTGCATCCACGCCCCGGTGACCTGCCGTTTGTTCAATAATTGACTGTGCCGGATCGCTGTCTTCATCAAAATTAATCGGGATCGCGCCGTAGCGGTCGGCGGCGAAATGCAAGCGGTAGGGATGATGATCAACAACAAAAATCTGTTCCGCACCGAGCAACCGTGCACAGGCGATTGTCAACAATCCCACAGGACCAGCACCATAGACTGC
AACGCTTGAACCTTGTTGGATCTGCGCATTTTTTGCTGCCTGCCATGCCGTTGGCAGAATATCAGAAAGGAAAAGCGCTTTATCATCTGAAAGCAAAGGCGGTACTTTAAACGGCCCCACATTCCCTTTAGGGACGCGGACATATTCCGCCTGCCCACCAGGAACGCCGCCATACAGGTGACTATAACCAAACAATGCCGC

which is the result of the assmebly using Megahit. While trying to read, it only recognizes k119_5 as the id, while I want the whole first line as the id of the sequence. Is vital information I want to keep!
Any ideas how to work with this?

The description field in the SeqRecord object has the information you are looking for:
>> from Bio import SeqIO
>> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file
>> s.id
k119_5
>> s.description
k119_5 flag=0 multi=141.0706 len=473

Edit: as an aside, if you write a SeqRecord object using SeqIO's write method, the FASTA header should contain the whole original header.


text:
True bam file encoding and viewing as text
We all have been using samtools to convert between sam and bam files.  SAM is tab-delimited file. My question resides on the ability to visualize bam files as text.  Is there any way to visualize any bam file as text in linux?
Why not work only with sam files and we need to convert them to bam?

What do you mean by "visualize bam files as text"? Do you search for something like IGV?

@Kamil S Jaron I was wondering how we can visualize as text, for example, less filename.sam would allow me to see SAM file as tab delimited.  Everyone has a good perspective on this topic.  For example @Pierre and @b.nota both mentioned that this is binary compressed file, therefore we can not see it directly without a program that interprets the bam binary encoding. An example of that program/interpreter is presented by @mgalardini and @gringer. And obviously, IGV being one of the widely used. Thank you, guys. Great answers.

Why not work only with sam files and we need to convert them to bam?

SAM is for human only, a binary format like BAM is smaller and parsed much faster by a program (and we want speed and optimize spaces for those huge files, don't we ?).
Furthermore, bam files are compressed using a  method that allow them to be accessed using random-access, that is how you can query a sub-part of bam file for a given genomic location.
you can always go back to SAM using samtools.

Assuming you want to visualize the pileup of reads across the reference, you may want to look into ASCIIgenome, a command line visualization tool for sam/bam/cram files, similar in spirit as IGV, but working in the command line.
ASCIIGenome my_sorted.bam

bam file is a compressed sam file, comparable with zip or tar. They are binary, which means we (humans) cannot read them. So you need to unpack them (convert back to sam) to read it as a text file.
bam files are smaller than sam files, which is the reason most people like to use bam files instead of sam files.

The way to visualise BAM files as text in Linux is with tools like samtools.
If you want to see the raw binary information as slightly more human-readable, you can run hexdump -C <file.bam>, but that type of text output is completely useless for further processing.
If you want more specifics about the BAM file format, you can have a look at the SAM file format specification, which also includes details on the related BAM and CRAM formats.


text:
Interpreting large z scores from microarray data
I have microarray z scores for gene set enriched (GSE) pathways that range approximately from -3 to 14. I understand that an absolute value of 1.96 has a p-value = .05 and significant for our purposes. My question is whether I can consider extreme z scores as more significant?

This is a statistical question. What a higher z value means is that it is more extreme (if you assume a normal distribution), thus a more extreme value is less likely to happen by chance. Which is translated to a (lower) p-value.
What it is significant and not is arbitrary. Prior to the experiment you must decide which is your threshold, how much chance you allow. Apparently you decided to select 0.05 as your threshold, which means that one of 20 results might be due to chance. But your extreme z score is not more significant, it is less expected (under the null hypothesis, those p-values might change if you change your null hypothesis, ie you think your genes expression don't follow a normal distribution) but as significant as your 0.0499999 p-values.


text:
Looking for datasets of RNA-Seq (or scRNA-Seq) for L1 larvae C. elegans
I was wondering if anyone know of a repository for C. elegans data sets, specifically in L1 larvae stage.
Are there any papers or GEO-sets for this organism at this stage?
thanks
Assa

See WormBase, which has a number of different resources for different nematodes:
https://wormbase.org


text:
Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?
For the purpose of testing phylogenetics software, we're using synthetic sequence data (i.e., it's computer generated, and there's no actual biological organism with that e.g. DNA sequence).  For real-world biological sequence data, we would say the sequence belongs to a "species" or "taxon" (or "taxonomic group"), but this feels inaccurate since there is no actual species or taxon.
Question: Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?
I'm wondering if this is standard practice, despite the synthetic nature of the data.

How were the synthetic reads generated exactly?

Using evolver, packaged with PAML.

Phylogenies can be made with/of non-species groups. I think "taxon" is perfectly fine. You might also consider the use of the "leaf / leaves" or "terminal" if you want to emphasise the tip nature of the organisms. Even "sequence" would be fine.
Ninja edit to my own suggestion: since you're making a phylogeny from sequences of unknown context or relationship (i.e. you can't say that they're sequences from a "species" or "OTU" or from anywhere), agnostic terms like "taxon" or "sequence" are more correct.

'Taxon' seems to be frequently used in the literature for simulated sequences. See this or this for example

When looking at genomic data for bacteria, taxa above the species level are operational taxonomic units (OTUs). OTUs are designated based on percent identity of the sequence (specifically the 16S rRNA gene sequences in most cases). The classification is usually based on rule of thumb, and some rough guidelines are:

97% for Genus
94% for Family
88% for Order

Using OTUs is common practice. I don't see much of an ideological difference using OTUs with synthetic data.
Konstantinidis, K. T., & Tiedje, J. M. (2005). Towards a Genome-Based Taxonomy for Prokaryotes. Journal of Bacteriology, 187(18), 6258–6264. http://doi.org/10.1128/JB.187.18.6258-6264.2005

OTUs are restricted to above the species level: they can be used for any group of closely related individuals. The use of OTU is meant to be agnostic as to actual group or relationship, underlining that the grouping was arrived at "operationally", i.e. methodologically. In any event, my feeling is that the context in which OTU is usually used is fairly distinct to a generic phylogenetic problem.


text:
Do additional peaks in percent GC of PacBio gDNA reads indicate contamination?
I have two sets of PacBio reads from genomic DNA of an Aspergillus species that were made from separate preps of  the culture. One of them has two additional peaks at 38% and 60% in the percent GC histogram produced by FastQC. Do these additional peaks indicate that there is contamination in that sample?

Yeah the prep looks odd, but you cannot really tell its contamination. You have to try and map the reads first, and see how many of them are unmapped, and see where the unmapped reads come from. You can just try with 1000 reads, should be obvious if it is contamination

Remake the graphs using only mapped reads. Then you can decide on your next steps.

The logic I would us is that if they come from the same genome, coverage peaks will go in multiplications of 2 - 1n, 2n, 4n, ..., while if it's a contamination, the coverage will go wild. So, I would recommend to look at the coverage vs GC to decide where it comes from.
You could either assemble the genome and use something like Blobtools to check coverage, GC and putative origin of sequences (and you should do that anyway).
Perhaps a bit simpler would be to look at it using k-mers, but that assumes you have reads with low error rate (In your case PacBio CCS or HiFi). If you do, you can use KAT that can do GC vs coverage plot - they intended it to look at GC bias in sequencing, but the vary same plot should show you if the three GC peaks have the same coverage or not.

I would suggest plugging your reads into a fast k-mer matching tool like mash screen to get a basic idea of where your reads come from, if you suspect contamination.
An issue with just looking at GC content is that some organelles/plasmids/etc. such as mitochondria have quite divergent GC content, and the different populations you observe might just indicate random fluctuations of copy number of these DNAs due to different sample preps. Not sure that would explain your plot but it does give me pause.
As suggested above, assembly is obviously a good idea if you are ok with the compute cost, and then you can use something like SPAdes or blobtools to pull out contigs from different origins.
KAT and similar tools as suggested are handy for hypothesis free evaluation of different populations in your reads, but they don't directly address the issue of contamination.
And obviously measure the quality of the reads to establish that they are ok and nothing weird happened there, they aren't all Gs or something.

Yes, it looks odd, but you might need to assemble the sequences to tell what is happening. The fact that the peaks are symmetrical around the presumably Aspergillus peak around 50% GC is strange. Usually contamination appears as one peak or shoulder on either side of the host DNA peak rather than symmetrical peaks.


text:
Using BLAT command line tool to blat split sequences
I have the nucleotide sequence:
AATTGAGGCACATTTTTTTTTAGACAGTCTTGCTCTGTTGCCCAGGCTGGAGTGCAGTGGTGTGATCATAGCTCACTGCAGCCTCGACCTCCTGGGCTCAACAAAGCACACAGTGGGCGGATCCCCACCAG
When I blat this on UCSC Genome Browser, the first hit is a match which spans 6572 basepairs on chromosome 19, and essentially matches the first portion of the sequence to one spot on the genome, and the second portion to another part of the genome approximately 6.5kb downstream.
However when I use the BLAT command line tool, it doesn't find the second half of this match, it only matches the first portion.
Here are the flags I'm using:
-t=dna
-q=dna
-minIdentity=90
-out=blast8
How do I get the command line tool to mimic the behavior of BLAT on UCSC Genome Browser?

The FAQ offers an answer:

I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.
We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. The web-based blat is tuned permissively with a minimum cut-off score of 20, which will display most of the alignments. Other than to confirm that your command-line blat is working, there is little use in perfectly replicating the web-based blat results. We advise deciding which filtering parameters make the most sense for the experiment or analysis. Often these settings will be different and more stringent than those of the web-based blat. With that in mind, use the following settings to replicate the search results of the web-based blat:
[ . . . ]
standalone blat:
blat search:
blat -stepSize=5 -repMatch=2253 -minScore=0 -minIdentity=0 database.2bit query.fa output.psl


text:
SnpEff not correctly annotating multiple-nucleotide polymorphisms
I have some SARS-CoV-2 sequencing data that I'm trying to annotate with SnpEff, however SnpEff doesn't appear to be recognizing multiple adjacent SNPs within the same codon and correctly calling them as a single MNP (multiple-nucleotide polymorphism).
I'm using the GenBank reference MN908947.3. Here are the relevant lines of the annotated VCF file
MN908947.3  28280   .   G   C   2911.06 PASS    AC=2;AF=1.00;AN=2;DP=65;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.56;SOR=8.380;ANN=C|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.7G>C|p.Asp3His|7/1260|7/1260|3/419||,C|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1278G>C|||||1278|,C|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2896G>C|||||2896|,C|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2060G>C|||||2060|,C|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1808G>C|||||1808|,C|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1089G>C|||||1089|,C|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*893G>C|||||893|,C|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*521G>C|||||521|,C|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*21G>C|||||21| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2925,196,0
MN908947.3  28281   .   A   T   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.91;SOR=8.380;ANN=T|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.8A>T|p.Asp3Val|8/1260|8/1260|3/419||,T|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1277A>T|||||1277|,T|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2897A>T|||||2897|,T|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2061A>T|||||2061|,T|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1809A>T|||||1809|,T|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1090A>T|||||1090|,T|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*894A>T|||||894|,T|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*522A>T|||||522|,T|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*22A>T|||||22| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0
MN908947.3  28282   .   T   A   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.71;SOR=8.380;ANN=A|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.9T>A|p.Asp3Glu|9/1260|9/1260|3/419||,A|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1276T>A|||||1276|,A|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2898T>A|||||2898|,A|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2062T>A|||||2062|,A|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1810T>A|||||1810|,A|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1091T>A|||||1091|,A|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*895T>A|||||895|,A|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*523T>A|||||523|,A|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*23T>A|||||23| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0

These are the positions corresponding to codon 3 on the N gene, they should be annotated as a single MNP, p.Asp3Leu, but istead they've been annotated seperately, as three individual SNPs.
I'm using the latest version of SnpEff (v5), and according to the manual SnpEff should have this functionality, what am I doing wrong?
You can see an example of how it's supposed to work here.

Are they expressed as a single variant in the input file? Does the manual say snpEFF can detect MNPs or only that it can handle them?

@terdon no, in the input file they're expressed as seperate SNPs. Based on the example in the link I provided, my take away from the SnpEff manual was that adjacent SNPs on the same codon would be grouped as an MNP, however there is little elaboration so I can't say for sure what is required for SnpEff to call adjacent SNPs as an MNP. I should note that since the variant caller does not have access to the ORFs, I don't think it would be able to call MNPs, so it has to be snpEff that does it.

The ORFs are not relevant here. In order to call a MNP, the variant caller needs to find all alt residues supported by the same read. Given that in your example, the reference allele has 0 support in all three, it seems very likely that the caller would find all three in the same read, but I don't see how snpEFF can be expected to change the list of variants you give it. If you join these yourself, snpEFF will likely deal with it, but I don't see anything in the manual that suggests it can merge them for you.

@terdon The ORFs are absolutely relevant here. The purpose of calling adjacent SNPs as a single MNP is so that the amino acid mutation is correctly annotated, thus you need to know where the breaks are between each codon, and cannot be deduced by the variant caller just from the ref. genome. Again, if you click on the link that I provided, and take a look at the example in the SnpEff manual, you will see how SnpEff annotates separate adjacent SNPs as a single MNP.

There is of course the problem of needing to determine phase (in the case diploid organisms), or even in the case of haploid organisms, that, as you say, SnpEff cannot in general determine from the VCF file whether adjacent SNPs always occur together on the same read. However I can find no info in the manual about how SnpEff handles these types of ambiguities, and my attempts at modifying the INFO field in the VCF file to indicate ploidy has not fixed the issue.

No, the ORFs are not relevant because the variant calling has no concept of genes, let alone translated products or codons. The calling is done using the reference genome and the reads. Of course the ORfs are relevant when you try to figure out the effect of your variants, but they do not affect and cannot affect whether something should be considered one variant or many. That is a function of the reads alone. All I see in the manual is an example of MNP, not a case where the input VCF had them separate and snpEFF joined them. I may be missing it, of course, but that's all I can find.

Let us continue this discussion in chat.


text:
Extracting specific sections from FASTQ determined by positions
Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions.
For example, say I want to extract the nucleotides at position 1, 3, 5-7 from each sequence:
Input:
@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

Output:
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)
Ideally, a command should work this way:
subselect 1,3,5-7 input.fq > output.fq

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5` and pasting) but maybe out of an box solution does already exist which evaded me.

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:
$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.
IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):
#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:
perl foo.pl file.fq 1,3,5-7

This looks like a job for samtools fqidx (http://www.htslib.org/doc/samtools-fqidx.html), assuming that your file is not enormous... read the warning in the link above.
usage would be
samtools fqidx yourfile.fq <region>
where region is specified using the fastq record identifier and coordinates.

How would a region like what the OP wants be written? You can run samtools fqidx file.fq id1 to get the data for id1, but how would you specify only the 1st, 3rd, 5th, 6th and 7th residue as requested by the OP? The docs don't show the format. And wouldn't this need each id in the fastq file to be mentioned individually? How would you use this for what the OP asked for?

I second terdons, fqidx is certainly a useful tool, but it doesn't address the specific question here.


text:
How are UMIs used to dedupulicate in Drop-seq tools?
The module DigitalExpression which is part of the popular Drop-seq tools digitally count gene transcripts.
The manual is not very clear on how exactly it resolves the UMI information.
Maybe I'm missing something but I was wondering if UMI are used in isolation or in conjunction with the barcode and or transcript?
The manual i.e. Drop-seq Alignment Cookbook only states:
To digitally count gene transcripts, a list of UMIs in each gene, within each cell, is assembled, and UMIs
within edit distance = 1 are merged together. The total number of unique UMI sequences is counted,
and this number is reported as the number of transcripts of that gene for a given cell.

To clarify my point, I'll use an example. Say my reads are only 4 bp and read 1 has the cell barcode at position 1+2 and UMI at position 3+4.  read 2 constitutes the transcript read.
Given the following example reads:
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = AATT  read2 = TTTT

UMI and barcode are the same but read2 is different, would they be collapsed to 1 UMI count or 2?
Similarly as case with same UMI and transcript but different cell barcode:
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = GGTT  read2 = ATGC

Also of note, UMI deduplication is the default setting as of the parameter below:
OUTPUT_READS_INSTEAD=Boolean  Output number of reads instead of number of unique molecular barcodes.  Default value:
                              false. This option can be set to 'null' to clear the default value. Possible values:
                              {true, false}

UMI and barcode are the same but read2 is different, would they be
  collapsed to 1 UMI count or 2?

If the two reads share cell barcodes and UMI sequences and align to the same gene, (even if they are different sequences) then they originate from the same RNA molecule, and should only count as 1.
https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry

Do you have any sources (documentation, code, citation etc) to back this up? To my knowledge the UMI is designed to spot PCR duplicates, so read 2 shouldn't have a different sequence (unless there are sequencing errors).

If you look up the library prep protocol, you can see that cDNA amplification happens before fragmentation.  Therefore, you can have multiple reads starting at different fragmentation points which nevertheless originate from the same mRNA molecule.  It's also quite plain from looking at the UMI counts output and bam output that all reads with the same gene ID, UMI and cell barcode are marked as duplicates save one, even if the reads aligned to different places in the transcript, and that only the one not marked as a duplicate is being counted in the final UMI count

Thanks. I Just had another look and this appears to be indeed the case, also for other drop-seq protocols apart from 10x. To however accept this as an answer, the only thing missing is if DropSeq-tools is doing the right thing by default or if this can be adjusted that way? The docu is surprisingly unclear about this important point.

It's been a while since I looked at DropSeq data, or its library prep protocol, but one has to presume that the Drop-Seq cookbook protocol handles DropSeq data correctly.  The description cited does seem to indicate that all the reads aligning to a gene are assessed for UMIs, not just all the reads aligning to a particular position.  However, I have in the past run the same Drop-Seq data set through the cookbook and (with some manipulation) through 10XGenomics cellranger, and I got very nearly the same results.

To confirm if this is correct, I've just submitted an issue

This is likely correct also for DropSeq tools and DropSeqPipe but this source is specifically for 10x Genomic Chromium platform. The source code for UMI processing for DropSeq tools can be found here: https://github.com/broadinstitute/Drop-seq

The code does EXACTLY what it says in the first quote.
Let's say you have 5 reads on a cell barcode that align to gene A.  Those reads have the following molecular barcodes:
AAAA,BBBB,AAAA,CCCC,DDDD.
There are 4 UNIQUE molecular barcodes.
The only complication to this is that we collapse barcodes with an edit distance of 1, which takes care of over-inflation of UMI counts due to PCR/Sequencing errors.  If we added a barcode with the sequence AAAB to the original list, it would be collapsed into AAAA as it has an edit distance of 1, and you'd still have 4 unique molecular barcodes.

In retrospect, the quote makes now sense, I was a bit confused what assembled exactly meant in this context, also the additional rational from swbarnes2 about the fragmentation taken place after PCR now explains why it is done this way (hence the accepted answer). Maybe this could be included in the documentation as explanation?


text:
How to segment genome into homozygote/heterozygote blocks based on RNA-Seq data
Given a RNA-Seq data of a progeny between two inbred lines (for example an F4 plant which was derived from two inbred lines (parents) and the F1 selfed 3 times; however this is probably not crucial for the question).
We know the progeny should have only a few breakpoints in each chromosome, so one should expect long stretches (>1MB) with those 3 states:

homozygous for parent1
homozygous for parent2
heterozygous

In each state, most containing genes should therefore be reflecting this state, i.e. only expressing the parent1 allele in state 1) (i.e all SNPs support parent1) etc.
A caveat would be uni-allelic expression at a heterozygous state 3), but this surely will not be the case for all genes in a given block (100+ genes per block are expected).
Both parental genomes and RNA-Seq libraries for both parents and F4 are available.
Naively I'd do SNP-calling (e.g. GATK HaplotypeCaller or bcftools call) of an F4 mapped against one parental reference and somehow combine the allele frequency in a sliding window and evaluate the likelihood of the given window to be in state 1-3. E.g. if all alleles support parent 1, it's probably state 1).
But I feel this must have been done before but can't quite find the right software or workflow software to use. But maybe I'm missing something.
Ideally, a workflow would have as input the F4 RNA-Seq data (possibly alos parental RNA-Seq data), the parental genome(s) and as output a list of segments per chromosome.
I know, we could do gDNA sequencing for this even better, but I think RNA-Seq might also do the job and might save effort since already available.


text:
What is the significance of the effect size on WOMAC scores for TKA?
https://www.medicaljournals.se/jrm/content/html/10.2340/16501977-2210
I understand the difference between statistical significance and effect size, but I am having trouble interpreting the data on Table III in the article linked above. I am evaluating the WOMAC scores. My first question is can effect size confidence intervals cross zero and still be valid? My second question is about the negative and positive values of the Hedges'g. The first study scores the WOMAC in a way where the higher the score, the better the outcome. The second two studies score the WOMAC in a way where the lower the score, the better the outcome. Does this mean that in the first study a positive Hedges'g indicate larger effect for the treatment group and a negative Hedges' g for the other two studies indicate larger effect for the treatment group?

Hi and welcome to the site. Please reproduce the relevant part of the paper here so we can understand what you are asking.

higher score = higher function Hedges'g g of -0.63 CI95% (-23.9, 5.9)

lower score = higher function. Hedges' g of 0.18 & CI95% (-2.9, 4.3)

Hedge's g measures the difference between two groups here experimental and  control in terms of combined (weighted) standard deviation. A higher g implies greater difference and lower g means lesser difference between  two groups.
It depends  on the study how the scoring of questions is proceeded with : the effect-size will be same. For example, answer to question item may be scored 1 to 10 depending on whether question item is written in positive  sense (the apple is good for health) or negative sense (the apple is not good for health). Positive Hedge's g does not indicate larger effect for the treatment group. THE negative Hedges' g for the other two studies indicates positive effect for the treatment group if the difference has been computed by  deducting experimental group mean from the control group mean (which is less than exp. Group Mean). Overall, I find that a knowledge of questionnaire scaling and scoring of scales and question items could be very helpful in general and in doing a meta-analysis effectively.


text:
Bacterial genome annotation of a clinical isolate strain?
So I'm basically so new to this that I'm just trying to find out what tools, methods, and keywords I should go look up by myself.
I have a unique strain of a bacteria.
I was given RNAseq data for this unique strain.
I want to analyze the RNAseq data, but I need an annotated genome to do that.
As a resource, I do have the annotated genome of the common lab strain.
The common strain should be pretty similar to the unique strain, but I don't know how similar.
Can anyone tell me what tools/methods might be applicable so I can go look them up?

For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain.
If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully de-novo, and one genome-guided. These transcripts can then be used as a reference assembly for Kallisto/Salmon, and counts compared to the lab-strain results. Note that the Trinity assemblies will only apply to the particular environments that are in your samples; if a transcript isn't expressed, it won't be in the assembly.
Any transcripts in the de-novo assembly that are not in the genome-guided assembly are potentially novel to your strain. However, care must be taken in interpreting this: if there is any other contamination in the sample, those contaminant transcripts will also be assembled.
I'm not too familiar with annotation pipelines, but based on a couple of Twitter posts it looks like Prokka might be a reasonable start. Torsten Seeman (lead developer of Prokka) has done a post about alternative annotation pipelines:
http://thegenomefactory.blogspot.co.nz/2013/03/bacterial-genome-annotation-systems.html
Update: Torsten Seeman now recommends Bakta as a "worthy successor to Prokka":
https://github.com/oschwengers/bakta

Bakta is a tool for the rapid & standardized annotation of bacterial genomes and plasmids from both isolates and MAGs. It provides dbxref-rich, sORF-including and taxon-independent annotations in machine-readable JSON & bioinformatics standard file formats for automated downstream analysis.

bacterial person here: I can confirm that prokka is the de-facto standard annotation pipeline for bacterial genomes, and should do a great job in annotating all gene features

Fantastic! Thank you so much for the great and detailed post, gringer! That helps so much! And thanks also mgalardini for confirming the info given!


text:
Bioconductor, how to retrieve the genes from a GEO gene expression dataset?
I'm quite new to Bioconductor and GEO.
Here's my task: I've to develop a simple R script able to read a gene expression dataset on GEO and extract the genes related to it. But... I don't know the command to extract the gene list.
The GEO dataset I want to analyze is GSE59867.
Here's my R code:
options(stringsAsFactors = FALSE)

library(oligo)
library(GEOquery)
source("https://bioconductor.org/biocLite.R")
biocLite("affyio")
library("affyio")

getGEOSuppFiles("GSE59867")
list.files("GSE59867")

untar("GSE59867/GSE59867_RAW.tar", exdir = "GSE59867/CEL")
list.files("GSE59867/CEL")

celfiles <- list.files("GSE59867", full = TRUE)
rawData <- read.celfiles(celfiles)

Now I have the rawData containing the gene expressions.
How can I extract the list of genes associated to rawData at this point?
Thanks!

Have a look at the pd.hugene.1.0.st.v1 package. That contains the annotations for that microarray.


text:
What is the R/Bioconductor preferred way of accessing TCGA data?
It seems there are two packages providing access to TCGA data these days:

TCGAbiolinks (older)
GenomicDataCommons (newer)

Will GenomicDataCommons supersede TCGAbiolinks? Is there anything that I might miss if I decide to use only GenomicDataCommons?

There isn't a R/Bioconductor person or view, so there isn't a preferred way to access TCGA. Use what suits you.
Bioconductor is build by contributors which are encourage to not overlap in their contributions. Usually the packages have different options/features even if they allow to access to the same data.
In this case, the TCGAbiolinks package is specifically for TCGA dataset while GenomicDataCommons is for "access the NIH / NCI Genomic Data Commons RESTful service". So the functionality they bring to R and Bioconductor is different.
The fact that one project is finishing doesn't mean that the data or the packages build around that project are obsolete. Some functionality different between them is the functions to harmonize and prepare (GDCprepare) data from the TCGAbiolink package, or see a boxplot correlation and AAIC plot to define outliers TCGAanalyze_Preprocessing.

I understand that there is not a R/BioC person, but it seems that the TCGA project is coming to a close in 2017. On the other hand the NCI GDC is just starting and will include/already includes the TCGA data, hence my question... Is it really true that they bring different functionality? Because it seems the GenomicDataCommons package will include the same (or more) of the TCGAbiolinks functionality.

@plant See the updates in the answer


text:
Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type?
Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type? For example, a FASTA file has 2 sequences. Can the first one encode amino acids while the second one encodes bases?
Thank you

While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes).
A better question to ask would be "Does ultra-cool-bioinformatics-tool X support combined nucleotide and protein sequences in the same FASTA file?" In which case the answer would most likely be, "No."

Yep. The FASTA format is very forgiving. Tools that consume data in FASTA format are...less forgiving.

If, on the other hand, Your software is dumb and will try to read aa sequence as nucletides, "fun" things will ensue.
Might I suggest using nexus format in this case?

The only obvious use case I can think of is simply using the fasta file to store all of your sequences and then extracting whichever you want to process from it. Any tool that can extract specific sequences from a multi fasta file shouldn't care whether those sequences are protein or nucleotide.

@terdon technically yes, but it's very bad practice. You literally wouldn't be able to hand this file over to someone else without additional metadata. There is no standard to specify whether the sequence is DNA or protein, and therefore no software package should be expected to be able to tell the two apart.


text:
What are the algorithms for biometric identification in large database?
I have been thinking about this problem and reading some related papers for sometime.
The problem can be described as follows:
Suppose that each biometric template is a $d-$dimensional vector. A database is composed of $n$ such templates ( $n$ could be large, like 1 million). Given a query template $q$, the task is to identify the "closest" template according to some distance metric.
Now, since the dimension $d$ and the database size $n$ could be both large, it is not quite efficient to linearly scan the whole database.
So my question is, are there any algorithms to efficiently fulfill the above task? What are the state of the art algorithms for biometric identification?
Through some literature review, I can think of two approaches. But both of them seem to have some limitations.

Building space partitioning index like R-tree. The problem is that such index does not suit for high dimensional data. The final search efficiency would be close to linear search.
Locality sensitive hash. The problem is that this method can only provide approximate searching result. For biometic identification, it is not guaranteed that the matched template could be found.

If this is properly stored in a database using indexes might be quite fast to retrieve the $n_i$ template if one variable should be identical. Otherwise it can be reduced to a n-dimensional through PCA techniques, and predict where the query q might fall on the dimensions, thus look on the nearest templates.

In practice: a) have cluster precompute all pairwise distances, so that you can query similarity to other records b) exploit specifics of your distance metric [with a million records or so this is usually still fast enough if your workstation is large enough to hold full matrix in memory, so that speed doesn't become an annoyance]; note that the possible modifications on the space and strategy and valid simplifications will depend on your distance metric

I think that a very nice example of what you are talking about saying the word bigger, is probably a microarray transciptional analysis of RNA extracted from a biological sample. Simplifying, in these analyis, each observation is analised on, up to 10000 genes expression levels of different mRNA coding genes. Thus, from these type of molecular biology experiments, we obtain matrices with few observations (d, the samples) respect to the dimensions observed (k, the target gene expression levels). Well, d << k. In this situation, in order to perform a classification-like task, and to reduce the problem of the overfitting related with dimensionality maledition, the ways that we might follow are several:
1) PCA: it is a good tool to visualize multidimensional datasets reducing dimensionality and noise level in the consolidated and standardized data. Problem: it is an unsupervised machine learning algorithm. Alone, it haven't the capacity of class discrimination.
2) PCA-LDA is a good combination, where PCA eliminates  hyperdimensionality and background in our data sheets, and than LDA (linear discriminant analysis) finds the directions onto the new subspace of characteristics, that are able to maximize variances between classes and decrease the variances within these. For this reason LDA is known as a kind of supervised machine learning technique. Furthermore, PCA and LDA are usually used together because LDA applied on high dimensional datasets could lead us toward overfitting problems. After the training of PCA-LDA pipeline, being LDA a classificator object trained and shaped on the principal component subspace, you can use this to examine new unknown samples or observation and predict its class by similarity. e.g. recently, Martin et al., on a Nature article talk about PCA-LDA approach to identify and than to predict cancer cells from those are normal, starting from ranman spectra of tissues.
3) PCA-LDA on transpose matrix: to use as better as possible the huge dimensionality, some authors used to transpose the initial matrix in order to force positioning of genes as rows and samples as columns. Then, the subsequent LDA is performed not on PC scores, but on loading factors of the eigenfactor matrix obtaining thus the list of the most important discriminant gene expression levels to segregate normal samples to cancer ones.


text:
snakemake rule all
I am trying snakemake for my mRNA-seq pipeline and I got stuck in a simple thing.
from os.path import join

RSEM_INDEX = "/app/ref/ensembl/human/idx/GRCh37"
SAMPLES, = glob_wildcards("Align/{sample}/Aligned.toTranscriptome.out.bam")

print("---Sample Names---")

for i in SAMPLES:
    print(i)

rule all:
    input:
        expand("Expression/{sample}/rsem.genes.results", sample=SAMPLES)

rule expression:
    input:
        bam = join("Align", "{wildcards.sample}}", "Aligned.toTranscriptome.out.bam")
    params:
        rsem_index = RSEM_INDEX
    output:
        join("Expression", "{wildcards.sample}", "rsem.genes.results")
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end Align/{sample}/Aligned.toTranscriptome.out.bam {params.rsem_index} Expression/{wildcards.sample}/rsem.genes.results >& Expression/{wildcards.sample}/rsem.log
        """

Here is an error I got
    MissingInputException in line 12 of /app/src/mRNAseq/TBD170049_20170320_RNAseq/quantify.snakemake:
Missing input files for rule all:
Expression/013_clamp_0m_buffy/rsem.genes.results
Expression/12_27_11_clamp_buffy_200/rsem.genes.results
Expression/010_200_buffy/rsem.genes.results
Expression/010_0_buffy/rsem.genes.results
Expression/12_27_11_clamp-p_buffy_0/rsem.genes.results
Expression/013_clamp_200m_buffy/rsem.genes.results
Expression/11_15_006_0m_buffy/rsem.genes.results
Expression/11_15_006_200m_buffy/rsem.genes.results

So, I thought rule all input for your final files you want to make, but it's looking for the final files that I want to make.
I would like to get some help. would you help me please? Thank you in advance!
My file trees.
drwxrwxr-x. 10 corea corea  4096 Sep 26 20:54 Align
-rw-rw-r--.  1 corea corea  1295 Sep 27 10:19 align.snakemake
-rw-rw-r--.  1 corea corea  1823 Sep 27 09:38 align.snakemake.test
-rw-rw-r--.  1 corea corea 31055 Sep 26 08:55 out
-rw-rw-r--.  1 corea corea   869 Sep 27 12:34 quantify.snakemake
drwxrwxr-x.  2 corea corea  4096 Sep 22 20:50 Rawdata
(mRNAseq) [corea@Biocore76 TBD170049_20170320_RNAseq]$ tree Align/
Align/
├── 010_0_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 010_200_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 013_clamp_0m_buffy

cross posted: https://www.biostars.org/p/274858/ and https://stackoverflow.com/questions/46439175

I am sorry. I will keep that rule in mind and it will never happen!

@jk_kim I don't think cross-posting is wrong per se, provided this is clearly acknowledged and that links are given so that readers of one site know there may be interesting information on the other site.

The following (or something close) should work:
rule expression:
    input:
        "Align/{sample}/Aligned.toTranscriptome.out.bam"
    params:
        rsem_index = RSEM_INDEX
    output:
        "Expression/{sample}/rsem.genes.results"
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end {input} {params.rsem_index} {output} >& Expression/{wildcards.sample}/rsem.log
        """

BTW, the original issue was probably due to the }} in your input line.

Thank you for the comment! no more cross posting! I will keep that rule in mind.


text:
Waterman-Smith-Beyer Sequence Alignment
I am having trouble understanding the affine gap penalty in the following example -

I am not sure where the 3 and 4 come from or the 4 and 5 in cells 1,3 and 2,4. I'm not understanding how the affine gap penalty works here. I understand how the linear gap penalty worked in the Needlemen_Wunsch but not this.

The linear gap penalty has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base.
The simplest case is assigning each unit (= each missing base) a penalty of 1:

g(k) = k

Or more generally, if you want to increase the penalty per base to 2, 3 or n:

g(k) = n * k
  with k = length of the gap, and n = cost of 1 missing base

The affine gap penalty has a cost that follows -as its name conveniently states- an affine function.
The idea behind it, is that any gap in an alignment has a cost with two components:

the "gap opening" penalty: it is a cost payed only once per gap, basically just for existing.
the "gap extension" penalty: equivalent to the linear gap penalty, it costs 1 unit for each consecutive missing base.

The general formula for the cost is then:

g(k) = a + n*k
  with a= opening penalty, n= extension penalty, k= length of the gap

Compared to the linear gap penalty, this will give a small bonus to longer consecutive gaps compared to skipping one nucleotide, matching one, skipping one again etc (because in the case of longer gaps, you only pay the opening penalty once, so the global cost per base decreases).
Cells on your diagram:
From what I understand, you're working with an affine penalty of g(k) = 1 + 1*k. For the [3,4] you're "coming" from the score of 2 on the left (cell 1,2). You either:

match the A (row 1) with G (column 3) and get a penalty of 1 for substitution -> total of previous 2 + new 1 = 3
open a gap to skip the G (column 3) and get a penalty of 2 (1 for opening + 1 because gap length of 1) -> total of previous 2 + new 2 = 4.

So you choose the lesser penalty and go with 3, which is the number selected in the bigger box.
The explanation is the exact same one for the other cell you mentioned.

I believe I understand now. The reason why we can get 3 or 4 is because in the previous cell, the optimal score of 2 could have come from either the left or the diagonal. If it had come from the diagonal, we are starting a new gap and thus we get 2 + (1 + 1) = 4. However, if we come from the left, we are simply adding to the already existing gap so we get 2 + 1 = 3. Is this the correct way to think about it?

The formula your slide shows describes a dynamic programming with a general gap cost. Finding the optimal solution takes $O(n^3)$ time. No one is using the algorithm these days. This is the first time I see an attempt to fill a matrix with this algorithm. It is interesting, but I am not sure how it works. Although affine gap cost is a special case of general gap cost, the algorithm to find the best score is very different, so is the time complexity. With modern affine gap formulations, each cell holds three values:

The overall optimal score (or the optimal score if this cell is a match, which yields an equivalent but different formulation)
The optimal score when the cell ends at a deletion
The optimal score when the cell ends at an insertion

The values in each cell only depends on the up, left and up-left cells in the matrix. You fill the matrix from the top left corner to the bottom right. This is an $O(n^2)$ algorithm, because you traverse each cell once. For more details, you can read this note. I would recommend to learn affine gap cost the normal way.

which algorithms produce local alignments under $O(n^3)$? Are there any R implementations of such algorithms?


text:
Sibling vs Parent/Child relationship detection using whole-genome sequence data
I am analyzing a family-based whole-genome sequence (WGS) dataset.  My identity-by-descent(IBD) analysis shows they share about 50% of their DNA. I think my pedigree structure may have gotten mixed up though and I want to see if one individual is the mother or the sister of a male individual.  Is there anyway to do this either with Plink or KING?

With at least three samples, you can determine this by looking at the IBS0 output column (from either plink2 --make-king-table or KING): if the relationship in question has a very low IBS0 value compared to the other pairs, that's parent-child.  You can try to do the same thing with only two samples, but with no other IBS0 values to compare against it's easier to misjudge.


text:
How to prevent sklearn Imputer(missing_values="NaN", strategy="mean", axis=0).fit_transform(data) from removing columns with only NA in them
I am trying to test a preexisting python machine learning script with a subset of my genetic data.  One of the feature columns I am using happens to only have NA values in it.  I lose this column when I use the following code:
Imputer(missing_values='NaN', strategy='mean', axis=0).fit_transform(data)
This causes problems later on because my model expects the correct number of features, and I do not want to retrain my model at this time.  Is there any way to work with Imputer or fit_transform to prevent this column from getting dropped?

Join the "NA" containing column after the transformation.

I am in no way expert in machine-learning but this seems more like a programming question than a bioinformatician problem, (although you use genetic data). Consider "telling" the model that there is now one feature less. But if you have new data enough different from the original you need to retrain your model AFAIK...

We might conceivably be able to help if you gave us some details but without your code or a dataset to test it on, we can't. You are using various commands there which, presumably, are defined in your model. We can't help if you don't tell us how they work or what they are.

authors: أماره أبو صهيون