cannin/biostars_qa huggingface dataset:
Question: I have access to some indexed BAM files in S3. Using the [AWS CLI][1], it is fairly easy to download the entire BAM file. That was fine for our initial analysis, but for validation, we are interested in looking a small region of ~10kb in thousands of individuals.

The BAM files are indexed and S3 supports GET requests with [range headers][2] so this should be possible. Does anyone know of a tool that does this?

EDIT: [htslib 1.3][3] was recently released and supports random access to BAM files in s3.

[1]: http://docs.aws.amazon.com/cli/latest/reference/s3/cp.html
[2]: https://docs.aws.amazon.com/AmazonS3/latest/dev/GettingObjectsUsingAPIs.html
[3]: https://github.com/samtools/htslib/releases/tag/1.3

Answer: The upcoming samtools 1.3 release will support this. As well as the current way of accessing public buckets in donfreed's comment on another answer, samtools 1.3 will understand s3: pseudo-URLs like

s3://1000genomes/phase1/data/NA12878/exome_alignment/NA12878.blah.bam

For accessing private buckets, samtools will look for your AWS credentials in the usual configuration files and environment variables, or you can specify them on the command line as *s3://id:secret@bucket/*... though that's not particularly recommended.

This release will be fairly soon. In the meantime, you can try this out by building samtools with GitHub htslib's **libcurl** branch. The code in that branch only looks for credentials in `$AWS_ACCESS_KEY_ID` / `$AWS_SECRET_ACCESS_KEY` and as *id:secret* in the URL.
Question: I found some tips for speeding up `read.table()` [here][1], I wondered if anyone can suggest something for [read.fasta()][2], which is part of the seqinr library. Not that they're really comparable.

I have a ~6GB file.

[1]: http://www.biostat.jhsph.edu/~rpeng/docs/R-large-tables.html
[2]: http://www.inside-r.org/packages/cran/seqinr/docs/read.fasta

Answer: If you have a huge file (just reread your post, 6Gb), the best way is to use an index (this is how databases work...). Take a look at this link: http://stackoverflow.com/questions/23173215/how-to-subset-sequences-in-fasta-file-based-on-sequence-id-or-name

It depends on what you want to do; do you only need a subset of these sequences at any given time or do you require every one of them? You would need a way to load and store them in your RAM, rather than reading a 6Gb file each time you tinker your file. NB: you can increase R's memory to avoid errors.
Question: Hi,

When I am asked to mention for which I have recently used a programming language like Python in data analysis, I just can say that I used Python once to extract some genome coordinates from gene bank. But definitely this is not a wise answer and disappointing for a bioinformatician. How then a medium data analyst could use Python? I mean for which purpose as a good answer.

Thank you very much

Answer: If this is hypothetical (and not explicitly about how *you* have used Python) then you can say that Python can be used for pretty much everything, although it has it's strong points and weaker points.

You could use python to create a pipeline (snakemake), make plots (matplotlib, seaborn), interface with databases (sqlite), perform processing of fasta/fastq (biopython) and create websites (django).

But even more important than knowing what you can is knowing the limitations. While you have a lot of statistical modules in Python - usually R is the more appropriate language here, definitely for differential expression analysis. But there is nothing that you can do in R that you cannot do in Python. For making production ready code for which speed is an issue then Python is often not your best guess, rather something compiled such as C or Java, or more recent languages such as Go or Rust. It's often 1) having a language you are comfortable in and 2) using the right tool for the job at hand.

I really dislike R, but I do differential expression analysis in R.
Question: Hello,

I'm analyzing TCGA breast cancer data to classify the samples into their respective subtypes, ad then to check if the genes of our study have a subtype specific pattern of expression. To do this, I was suggested to use `genefu`. At the step of classifying the subtypes -

`PAM50Base <- molecular.subtyping(sbt.model = "pam50",data=data, annot=annot,do.mapping=F)`

I get an error -

`Error in intrinsic.cluster.predict(sbt.model = pam50.robust, data = data, : `

`no probe in common -> annot or mapping parameters are necessary for the mapping process!`

In the command, annot is the file used for annotation and is of the format -

>probe EntrezGene.ID Gene.ID Gene.Symbol

Data refers to the input file, which is of the format -

>Gene.Symbol Sample1 SAMPLE2 ... Sample 1092

Both files are tab delimited. I want to know if anyone has done this before, and if the file formatting is correct?

P.S. I have tried using the `Gene.Symbol` and `probe` in the data file, but both give the same error.

Edit: Should my data file also contain the `EntrezGene.ID` column?

Thank you.

Answer: Hi vinayjrao,

I think the problem is with the `data` matrix; `molecular.subtyping` fuction expects a matrix of samples(rows) x genes(cols). As I can see above your data matrix is genes x samples, right? Try transposing the matrix.

From the genefu vignette for do.mapping: `TRUE if the mapping through Entrez Gene ids must be performed (in case of ambiguities, the most variant probe is kept for each gene)`
Question: Hello,

In our group we are about to design an amplicon-based targeted panel for NGS, and we're aiming at around 150 bp for mean amplicon sizes. The wet lab people came back to me asking to calculate the optimal number of samples to be loaded in the instrument (MiSeq). And this is where I have my question: ideally I'd like to pick 2x150bp as read length, so I can have a greater output, however the mean amplicon size might be close to the read length, or in some cases, lower.

Hence my question: would it be "safe" to use a read length approximately equal to the amplicon length so that it doesn't complicate my analysis job afterwards (I'm thinking about adapter read through)? My previous experience always had a larger amplicon length (> 200 bp) so it can't really compare.

Should I still go for 2x150bp, or aim lower at 2x75bp? If I use 2x150bp, would possible adapter read through throw a spanner in the works, or should soft clipping (from BWA) be enough?

Answer: I'd still do 2x150 if the amplicon targets are 150bp on average. You can deal with the adapter read through and trimming fairly easily. Something like 200bp averaged amplicons might be a better size overall, but either way I think is fine.
Question: Trying to understand substitution matrices. It seems like it is a scoring scheme for alignments, particularly if you are looking for homology?

I am trying to see if it would be applicable if I am looking at mutations between proteins from different people. Since my sequences are very similar with only 1 or 2 mutations between them, the substitution matrix would probably not be applicable here? I am assuming if there is a nonsynonymous mutation between two sequences it would give me a score (say BLOSUM62) based on how likely that substitution would occur in nature? Are there other ways to interpret these scoring matrices?

Answer: Before we go into your question, it may be best and most concise to simply describe the exact SNP sites and leave it at that, given that your proteins are so similar. However, here are the differences in PAM and BLOSUM:

BLOSUM (BLOcks SUbstitution Matrix) were derived by looking at alignments of highly conserved protein domains at different evoluntionarily divergent distances, then taking into account how frequently one amino acid was substituted to another. It's described in [this paper by Henikoff][1]. They are based on local alignment of conserved protein regions.

PAM (Point Accepted Mutations) matrices were first described by Margaret Dayhoff (who was a fantastic scientist, even in face of the challenges of her role given the time period). "Each entry in a PAM matrix indicates the likelihood of the amino acid of that row being replaced with the amino acid of that column through a series of one or more point accepted mutations during a specified evolutionary interval, rather than these two amino acids being aligned due to chance." They are based on global alignment.

In short, this is what matters about the differences between the two:

1. PAM matrices are typically used on more closely related proteins (such as your case), BLOSUM are typically used on more evolutionarily divergent proteins.
2. The greater the PAM number the more DISTANT the sequences being compared should be; the greater the BLOSUM number, the more SIMILAR the sequences being compared should be.

So for your application, if you were to use these, you should either use a LOW PAM matrix or a HIGH BLOSUM matrix number. Whether this is appropriate for your application depends on what you want to get out of it (e.g. the whole protein difference or just local protein domain differences); you're right in that they are typically used for alignment scoring, but they can also be used to generate some evolutionary cost distance. However, there may be better methods out there for your purpose if you look for methods for creating distance trees based on some metric.

[1]: https://www.ncbi.nlm.nih.gov/pubmed/1438297
Question: Hi there,

I always used TCGAbiolinks to get raw count for TCGA projects like below:

expquery <- GDCquery(project = "TCGA-KIRC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
GDCdownload(expquery,directory = "GDCdata")
expquery2 <- GDCprepare(expquery,directory = "GDCdata",summarizedExperiment = T)
expMatrix <- TCGAanalyze_Preprocessing(expquery2)

However, it does not work today and it seems there is no HTSeq - Counts

Error in GDCquery(project = "TCGA-KIRC", data.category = "Transcriptome Profiling", :
Please set a valid workflow.type argument from the list below:
=> STAR - Counts

Therefore, I used STAR - Counts to download the data, which has completely different format of what I downloaded before using HTSeq - Count. The expression matrix for each sample has more columns including fpkm_unstranded and tpm_unstranded.

Actually the data downloaded using STAR - Counts is much more useful but I do not know how to extract the files to a readable expression matrix (ENSEMBL ID as rownames and TCGA tumor barcode as column names) because the GDCprepare() function fails to work on it:

expquery2 <- GDCprepare(expquery,directory = "GDCdata",summarizedExperiment = T)
| | 0%
Error in readr::read_tsv(file = f, col_names = TRUE, progress = FALSE, :
unused argument (show_col_types = FALSE)
Error in if (value == n) { : argument is of length zero

Anyone has any solutions? Many thanks in advance!

Answer: ```
projects <- c("TCGA-LGG","TCGA-GBM")

ids <- c("TCGA-76-6661", "TCGA-26-5132", "TCGA-19-1389")

# RNA Transcriptome
query <- GDCquery(
project = projects,
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts",
barcode = ids
)

GDCdownload(query,
method = "api",
directory = "./data",
files.per.chunk = 10)

rna <- GDCprepare(query, save = FALSE, summarizedExperiment = T, directory = "./data", remove.files.prepared = F)
rna <- as.data.frame(rna@assays@data@listData)

# DNA Meth Beta Value
query <-GDCquery(project = projects,
data.category= "DNA Methylation",
platform = "Illumina Human Methylation 450",
data.type = "Methylation Beta Value",
legacy = FALSE,
barcode = ids)

GDCdownload(query,
method = "api",
directory = "./data",
files.per.chunk = 10)

met <- GDCprepare(query, save = FALSE, summarizedExperiment = T, directory = "./data", remove.files.prepared = FALSE)

# IDAT Files
query <-GDCquery(project = projects,
data.category = "Raw microarray data",
data.type = "Raw intensities",
experimental.strategy = "Methylation array",
legacy = TRUE,
file.type = ".idat",
platform = "Illumina human methylation 450",
barcode = ids)

GDCdownload(query,
method = "api",
directory = "./data",
files.per.chunk = 10)
```
Question: I love .gfa, but sometimes I have trouble to understand them.

I have used Flye with pacBio reads with defaults options to make a first shot at assemble a linear bacterial genome. The genome probably contains some plasmids or some phages sequences. Flye gave me the following .gfa :

![enter image description here][1]

[1]: /media/images/4de3edca-7c14-4314-9dcc-e6f90b16

The green, yellow and red edges have been merged into one scaffold of around 9 Mb (as expected). And the little blue edge is his own contig of 42 Kb.

In his assembly_infot.txt file, Flye report that the big scaffold is indeed non linear, and that the tiny one is circular :

#seq_name	length	cov.	circ.	repeat	mult.	alt_group	graph_path
scaffold_2	8926751	84	N	N	1	*	*,1,2,4,??,4,-3,-1,*
contig_4	42310	940	N	Y	12	*	4

I have a few questions that puzzles me about this graph :

- Why does the .gfa connect all the edges into a circle if flye report that only one piece is circular not the other ?
- The yellow edge is connected by the same end to both green and red... And if the mean coverage of the green and red is 70X, it is only 17X for the yellow edge. What could it mean ? I am very puzzled by the fact it is connected by the same end. Could it be Flye trying to circularize it ? Or a sort of SV ? I think the DNA provided comes from a single colony so I don't see how that could be a SV.
- The blue edge has connection to itself, I imagine it is because of repetitions. But in the graph, this repetition is somehow connected to the others edges. So why was it split in the .fasta file in the end ?

Answer: Oh, I forgot to update this post, but fenderglass, the author of flye, helped me a lot to understand this. You can follow the whole conversation here : https://github.com/fenderglass/Flye/issues/389
Question: Hello Biostars,

What are an approximate quantity of human genes, transcripts (coding transcripts), proteins and metabolites? What about a model plant like Arabidopsis?

Thanks

Answer: Here are some statistics that I found in different databases:

For human:

- Coding genes: about 20,441
- Non-coding genes: about 22,219
- Transcription factors: about 2,067 (TFs/total genes= 10 %)
- Transcripts: about 198,002 (Transcripts/genes = 9.68)
- Proteins: more than 30,057
- Metabolites: more than 42162

Any revision or adding more information to this post would be
appreciated,
Question: I want to get the information of all genes on human Y chromosome, then I found the statistics in different databases --Ensembl (GENCODE), NCBI, HGNC -- are dissimilar.

For example, protein-coding genes numbers:

CCDS	63
HGNC	45
Ensembl	63
NCBI	73

So what leads to these number be different?

By the way, is RefSeq gene data the same as NCBI homo sapiens annotation release?

Answer: The different resources you cite do different things and are not necessarily in sync. The CCDS tries to identify annotations of protein-coding regions in the human and mouse genomes that are consensual across several groups/institutes. The HGNC is in charge of attributing official names and symbols to genes. NCBI's RefSeq is a collection of sequences that are annotated as belonging to a gene and/or linked to other NCBI resources.
Ensembl provides a full genome annotation integrating many information types. Of the resources you cite, Ensembl is the only one that annotates the underlying genome.
When doing a bioinformatics project, select one reference and stick to it. Don't mix and match, this would be asking for trouble. I would recommend using Ensembl because it's much better organized and integrated than NCBI resources.
Question: I'm analyzing gene expression data from an experiment where I have (i) 6 different mice (A:D); (ii) 2 different time points ("pre" and "post" treatment); and (iii) 2 different response types ("responder" vs. "non-responder"). I'm trying to find genes with differential expression after treatment in the responders vs. non-responders. In other words: if I look at post therapy gene expression changes compared to 'pre', which genes are differentially expressed in responders compared to non-responders? Additionally, how can I include the 'paired' aspect into my analysis, where I might get more meaningful information by running the above analysis in a paired, rather than general, manner?

I've spent an entire day reading tutorials on designing design matrices for differential expression but I can't pin down the details, I'm hoping someone with experience might be able to help out.

A small reproducible example:

ms = factor(rep(c("A", "B", "C", "D", "E", "F"), each = 2)) # Mouse ID
rx = factor(rep(c("Pre", "Post"), 6)) # Pre or Post treatment
rp = factor(c(rep("NR", 6), rep("R", 6))) # Is the mouse a 'responder' or 'non-responder'?
data.frame(ms, rx, rp)
ms rx rp
1 A Pre NR
2 A Post NR
3 B Pre NR
4 B Post NR
5 C Pre NR
6 C Post NR
7 D Pre R
8 D Post R
9 E Pre R
10 E Post R
11 F Pre R
12 F Post R

The ideal analysis that I'd like to do is something like the following, but including the paired design, which would hopefully make the analysis more powerful:

(Post / Pre in responders) / (Post / Pre in non-responders)

The first term ("Post / Pre in responders") would find DE genes in responders, and by dividing by the second term ("Post / Pre in non-responders") I would ideally get genes which are diferentially expressed only in responders.

I tried running the "Post / Pre" analysis in responders and in non-responders separately, and manually dividing, as in the code above, but it doesn't feel right, and there must be a true way of running that with a properly designed matrix.

Answer: You'll want a design of `~ms + rx*rp` and the coefficient of primary interest is then the `rx:rp` interaction. That then accounts for the pairing and gives you the interaction you want. Note that you'll probably get an error about the matrix being of insufficient rank, in which case you should use `ms = factor(c(rep(rep(c("A","B","C"), each=2), 2))`. This will keep the pairing, but change things such that each of the groups contains a mouse with the same ID (otherwise correcting for the pairing will confound fitting the response group).
Question: Hi!

I'm trying to run fastQC on an HPC system but i'm having this error:

It seems that something is wrong with java.

Started analysis of raw_f.fastq

Approx 5% complete for raw_f.fastq

Approx 10% complete for raw_f.fastq

Approx 15% complete for raw_f.fastq

Approx 20% complete for raw_f.fastq

Approx 25% complete for raw_f.fastq

Approx 30% complete for raw_f.fastq

Approx 35% complete for raw_f.fastq

Approx 40% complete for raw_f.fastq

Approx 45% complete for raw_f.fastq

Approx 50% complete for raw_f.fastq

Approx 55% complete for raw_f.fastq

Approx 60% complete for raw_f.fastq

Approx 65% complete for raw_f.fastq

Approx 70% complete for raw_f.fastq

ecc...

javax.imageio.IIOException: Can't create cache file!
at javax.imageio.ImageIO.createImageInputStream(ImageIO.java:361)
at javax.imageio.ImageIO.read(ImageIO.java:1397)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.base64ForIcon(HTMLReportArchive.java:379)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.startDocument(HTMLReportArchive.java:303)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.(HTMLReportArchive.java:84)
at uk.ac.babraham.FastQC.Analysis.OfflineRunner.analysisComplete(OfflineRunner.java:185)
at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:123)
at java.lang.Thread.run(Thread.java:745)
Caused by: java.nio.file.AccessDeniedException: /tmp/imageio6327552574161602354.tmp
at sun.nio.fs.UnixException.translateToIOException(UnixException.java:84)
at sun.nio.fs.UnixException.rethrowAsIOException(UnixException.java:102)
at sun.nio.fs.UnixException.rethrowAsIOException(UnixException.java:107)
at sun.nio.fs.UnixFileSystemProvider.newByteChannel(UnixFileSystemProvider.java:214)
at java.nio.file.Files.newByteChannel(Files.java:361)
at java.nio.file.Files.createFile(Files.java:632)
at java.nio.file.TempFileHelper.create(TempFileHelper.java:138)
at java.nio.file.TempFileHelper.createTempFile(TempFileHelper.java:161)
at java.nio.file.Files.createTempFile(Files.java:897)
at javax.imageio.stream.FileCacheImageInputStream.(FileCacheImageInputStream.java:102)
at com.sun.imageio.spi.InputStreamImageInputStreamSpi.createInputStreamInstance(InputStreamImageInputStreamSpi.java:69)

I installed the latest fastQC version inside the HPC on my personal conda environment but I'm still getting the same error.
Fastqc creates the report.txt file but the HTML file is missed and i can't see any plot in the zipped file generated.

Any suggestion?

Answer: your problem is here:

> Caused by: java.nio.file.AccessDeniedException: /tmp/imageio6327552574161602354.tmp

use option -d

-d --dir Selects a directory to be used for temporary files written when
generating report images. Defaults to system temp directory if
not specified.
Question: When analyzing clinical genomes, I often find it hard to visualize the genomic context of a potential driver mutation (say a point mutation on a key cancer gene). For example, is there simultaneous LOH/amplification at the mutation locus? How many copies of the mutant allele are there? Did the mutation arrive before/after the overlapping copy number change? Is there a structural variation breakpoint near by? Are there other point mutations close by?

Thus, I wonder if there is a plotting tool that visualizes a local range of the genome (e.g. a window big enough to inspect a gene) that integrates (in one plot):

1. Point mutations in the region
2. The local copy numbers (major/minor allele if available)
3. SVs affecting the region

Note that the purpose is more to infer the consequence of a mutation as opposed to validating whether the mutation is real or not, which we can already do with a variety tools such as samtools tview.

Circos is a nice format that integrates all three classes of mutations, but it is more for global scale inspections and does not provide enough resolution to zoom into specific gene regions in order to answer the above questions about key mutations of interest. Another close one is IGV, but it does not incorporate allelic specific copy number info.

Thanks!

Answer: You can give a try to [karyoploteR](http://bioconductor.org/packages/karyoploteR/). It's an R package to plot data on genomes and can go from single base to the whole genome. Take a look at the [tutorial and examples](https://bernatgel.github.io/karyoploter_tutorial/) to see what it looks like. It's only a plotting tool and does not know much biology, so it can plot almost anything as long as you have it precomputed and can map it to any of the available plot types. Right now it does not have the possibility of plotting individual sequencing reads from a BAM file, if you need that.

A couple of examples of whole genome views.

Note: Setting `zoom` to a small genomic region in the call to `plotKaryotype` will created a zoomed in version of the plots with the exact same code.

![enter image description here][1]
![enter image description here][2]
![enter image description here][3]

[1]: https://bernatgel.github.io/karyoploter_tutorial/Examples/SNPArray/images/Figure4-1.png
[2]: https://bernatgel.github.io/karyoploter_tutorial/Examples/GeneExpression/images/Figure13-1.png
[3]: https://bernatgel.github.io/karyoploter_tutorial/Examples/Rainfall/images/Figure5-1.png