cannin/biostars_qa huggingface dataset:
Question: Hello Biostars Community,

I am a newbie to the bioinformatics world and need some help with a project I am working on!

I need to use htsjdk to extract every tlen value from a BAM file and have a few questions.

- How do a read the BAM file in the first place? I have found SamReader that might do the job like so :
SamReader samReader= SamReaderFactory.makeDefault().open(new File("mybam.bam"));
- How do I find the position of individual chromosomes? I want to multithread the application, reading each chromosome in parallel.
- Once I get the position to start reading from, how do I read each line, then the tlen field?

Sorry for all of the questions, I have been googling a lot and haven't found answers!

Any help is appreciated, thanks in advance!

Answer: This may get you in the right direction:

// Open bam file
SamReader samReader= SamReaderFactory.
makeDefault().
validationStringency(ValidationStringency.SILENT).
open(new File("aln.bam"));

// Start iterating from start to end of chr7.
SAMRecordIterator iter = samReader.query("chr7", 0, 0, false);

while(iter.hasNext()){
// Iterate thorough each record and extract fragment size
SAMRecord rec= iter.next();
int tlen= rec.getInferredInsertSize();
System.out.println(tlen);
}

Get list of chromosomes/contig names from `SamReader` object:

List refSeqs= samReader.getFileHeader().getSequenceDictionary().getSequences();
for(SAMSequenceRecord ref : refSeqs){
System.out.println(ref.getSequenceName());
}
Question: Hi,

I have my bam file in this order:

samtools view -H my.bam

@CO Duplicates marked (bamUtil Version: 1.0.2)
@SQ SN:chr1 LN:249250621 AS:hg19
@SQ SN:chr10 LN:135534747 AS:hg19
@SQ SN:chr11 LN:135006516 AS:hg19
@SQ SN:chr12 LN:133851895 AS:hg19
@SQ SN:chr13 LN:115169878 AS:hg19
@SQ SN:chr14 LN:107349540 AS:hg19
@SQ SN:chr15 LN:102531392 AS:hg19
@SQ SN:chr16 LN:90354753 AS:hg19
@SQ SN:chr17 LN:81195210 AS:hg19
@SQ SN:chr18 LN:78077248 AS:hg19
@SQ SN:chr19 LN:59128983 AS:hg19
@SQ SN:chr2 LN:243199373 AS:hg19
@SQ SN:chr20 LN:63025520 AS:hg19
@SQ SN:chr21 LN:48129895 AS:hg19
@SQ SN:chr22 LN:51304566 AS:hg19
@SQ SN:chr3 LN:198022430 AS:hg19
@SQ SN:chr4 LN:191154276 AS:hg19
@SQ SN:chr5 LN:180915260 AS:hg19
@SQ SN:chr6 LN:171115067 AS:hg19
@SQ SN:chr7 LN:159138663 AS:hg19
@SQ SN:chr8 LN:146364022 AS:hg19
@SQ SN:chr9 LN:141213431 AS:hg19
@SQ SN:chrM LN:16571 AS:hg19
@SQ SN:chrX LN:155270560 AS:hg19
@SQ SN:chrY LN:59373566 AS:hg19
@PG ID:bwa PN:bwa VN:0.6.1-r104

And I wanna order it as, this:

https://gist.github.com/RamRS/510a730a4b77f8ded1d95acb97067e3c

How can I sort it in this order and arrange it?

Thank you

Answer: use [picard **ReorderSam**](http://picard.sourceforge.net/command-line-overview.shtml)

> ReorderSam reorders reads in a SAM/BAM file to match the contig ordering in a provided reference file, as determined by exact name matching of contigs
Question: Hello everyone,

I've tried to work with the TCGA DNA methylation data, but I'm having problems to understand better these data.

The TCGA's website isn't working anymore to get the biological data. In its place this GDC website (https://gdc-portal.nci.nih.gov/) is working. I was able to get clinical, mRNA and miRNA data in that site, however, I didn't find the DNA methylation data. Isn't there dna methylation data in that portal?

Fortunately, I found another site: cancer genomics browser (https://genome-cancer.ucsc.edu/proj/site/hgHeatmap/) where I was able to get the DNA methylation data for breast cancer (HumanMethylation450).

There are multiple files in the dna methylation file. levels of methylation are in the "genomicmatrix.txt" file, which each sample methylation has beta-value and a probe. On the other hand, the probe.txt file contains correspondence between the probe with the genes. Here is a little example about the genomicMatrixand probe file.

genomicMatrix:

sample	TCGA-OL-A66H-01	TCGA-3C-AALK-01	TCGA-AC-A5EH-01
cg13332474	-0.4808	-0.2968	-0.1997
cg00651829	-0.4821	-0.2110	-0.4108
cg17027195	-0.4633	-0.4250	-0.4667
cg09868354	-0.4345	-0.3630	-0.4230
cg03050183	-0.4252	-0.3749	0.1269

cg01989731	NA	NA	NA
cg06819656	0.4028	0.3047	0.3755
cg04244851	0.4398	0.3894	0.2533

cg19669385	-0.1353	0.3650	0.0664

cg04244855	0.4292	0.4008	0.2468

cg17689707	-0.4842	0.0109	-0.2484
cg04244857	-0.0918	0.2731	-0.0084
cg02434381	-0.4443	-0.4273	-0.4175
cg05777492	-0.4595	-0.4780	-0.4786
cg23340034	0.0933	0.3611	0.4120

cg26361545	0.4339	0.4389	0.4348

cg10609310	0.2913	0.0337	-0.1307

When looking at the file genomicmatrix.txt see several negative and NA values. I thought of disregarding them. Have positive values, I do not find any value above 0.8, ie no hypermethylation values. Why?

Probe:

id	gene	chrom	chromStart	chromEnd	strand

cg00035864	TTTY18	chrY	8613009	8613010	.

cg13275322	WAS	chrX	48426764	48426765	.

cg13798679	chr1	36390157	36390158	.

cg13799227	chr1	226719204	226719205	.

cg13799302	CYP2J2	chr1	60164980	60164981	.

cg13799671	CD58	chr1	116881090	116881091	.

cg13805052	MORN1,LOC100129534	chr1	2272923	2272924	.

Here I consider only genes that are not on chromosome X and Y. I noticed that there are some probes that associate with more than one gene, in this case thought to obtain the median of methylation values ​​to result in final gene methylation level.

I was thinking of converting these files into a single file with the following header:

gene | beta value | sample
CD58| 0.4 | TCGA-OL-A66H-01

Please can someone help me these questions?

Thank you for attention.

Answer: Hi!

Thank you very much for helping me!

I got the data in the GDC legacy web site. Unfortunately, there aren't many patients in both conditions, e.g. Tumor and Normal. I suppose that there are 40 paired samples and about 300 samples.

About the DNA methylation values, it makes sense the values are between -0.5 and 0.5, but can I apply the rule of three - proportions in these values? Or there is another manner to convert these values to beta values?
Question: I am confuse if padj and FDR is same? I am in the mid of analysis and confused if I have FDR or not. I have seen some explanation but not able to get the answer. please share your views.
Thanks

Answer: In DESeq2 the FDR-corrected p-values are called `padj` indeed.
Question: I have a fasta file which contains thousands of sequences, with headers as such:

>scaffold_1|c(135298..135582)|DNA|DNA-0-1_NV

Each pipe-deliminated section of the header can vary from sequence to sequence, and some sequences might have identical headers except for the first or second sections.

I need to be able to search through this large file and pick out and print to another file specific sequences based upon their header. There needs to be degeneracy in this search however. I have seen examples where a library text file was used but only exact matches between the fasta file and library file would work.

For instance, let's say I want all sequences which have any variation on 'piggyBac' in their header (so PiggyBac, piggybac,DNA-piggyBac, etc.).

I'm just at a loss as to how to do this exactly. Is there some way to index this file and then search the keys for variations on 'piggyBac'? If anyone has suggestions or can point me to code that does something similar it would really be helpful.

I appreciate it

Answer: The most elegant solution would use an iterator approach. This is based on the "Filtering a sequence file" example in [the tutorial][1] only rather than filtering based on the identifiers, here we are checking the description:

```
from Bio import SeqIO
input_file = "big_file.fasta"
output_file = "piggybac.fasta"
# The next line is a Python generator expression - memory efficient!
records = (r for r in SeqIO.parse(input_file, "fasta") if "PIGGYBAC" in r.id.upper())
count = SeqIO.write(records, output_file, "fasta")
print("Saved %i records from %s to %s" % (count, input_file, output_file))
```

Updated with variable name correction pointed out by [tyleraelliott][2]

[1]: http://biopython.org/DIST/docs/tutorial/Tutorial.html
[2]: https://www.biostars.org/u/15855/
Question: Hi:
I have a BED file with sites of a transcription factor binding sites (mm9). I want to annotate these sites, typically located between intergenic regions, to a gene nearby. I checked a publication (PMC3080774) and it lists some tools that can map and annotate.
I tried using Chipseqanno. Initially it was straightforward till first two commands (https://bioconductor.org/packages/release/bioc/vignettes/ChIPpeakAnno/inst/doc/pipeline.html) and it went confusing.

bed <- system.file("extdata", "MACS_output.bed", package="ChIPpeakAnno")
gr1 <- toGRanges(bed, format="BED", header=FALSE)
## one can also try import from rtracklayer
gff <- system.file("extdata", "GFF_peaks.gff", package="ChIPpeakAnno")
gr2 <- toGRanges(gff, format="GFF", header=FALSE, skip=3)
## must keep the class exactly same as gr1$score, i.e., numeric.
gr2$score <- as.numeric(gr2$score)
ol <- findOverlapsOfPeaks(gr1, gr2)
## add metadata (mean of score) to the overlapping peaks

Authors did not explain why they are overlapping gr1 and gr2. Before that authors say "## one can also try import from rtracklayer". also authors did not point out names for annotations MM8 or Hg19.

My question to this forum is:
1. Given a bed file how can I map and annotate the peak regions to nearby genes.
2. Can I define 1kb, 2kb range. (apparently chipseqanno uses granges and thus is possible- however lack of mm9 and proper reference to why findOverlaps between gr1 and gr2 makes it so confusing to deal with ranges).
3. Are there any tools other than chipseqanno to achieve above two tasks.

Thanks
Adrian

Answer: You can try using [annotatePeaks.pl][1] which is a simple perl script, a part of [HOMER][2] package. You basically need 3 files for running it:

1. Peak bed file
2. reference genome fasta file
3. reference genome annotation (gtf/gff) file

And command goes like this:

perl annotatePeaks.pl sample_peak.bed Genome.fasta -gft genome.gft >homer_annotation.txt

[1]: http://homer.ucsd.edu/homer/ngs/annotation.html
[2]: http://homer.ucsd.edu/homer/
Question: I am looking for a tool that can do a insillico digestion of a reference genome given 2 restriction enzymes (or the pattern at which both cut).

The output that I would like to have is a BED file with the start and end position on the reference genome of all the produced RE DNA fragments. (or only those that are larger than a specific lenght).

The tool should take into account that both the forward and reverse DNA strand could be cut.

Is there any such tool that I can download?

Answer: I found the R package SimRad which

> provides a number functions to simulate restriction enzyme digestion,
> library construction and fragments size selection.

https://cran.r-project.org/web/packages/SimRAD/SimRAD.pdf

A nice feature is that a filter step can also be done on restriction site combination to only select fragments that start with enzyme A restriction site and end with enzyme B restriction site.

adapt.select(sequences, type = "AB+BA", cut_site_5prime1, cut_site_3prime1,
cut_site_5prime2, cut_site_3prime2)

The output is the DNA fragments produced by the digestion and size selection.

The output does not contain were the fragments are located on the genome.

To create the BED file I need to BLAST the fragments versus the reference genome I guess.
Question: I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

```
>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence
```

And I have a list like the following:

```
Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11
```

I want to retrieve those sequences that match with names present in the list.

Answer: You could use Biopython's fasta parser. I didn't test this, but this would be the basic idea:

```py
from Bio import SeqIO
handle = open("example.fasta", "rU")
organism_name = "Paenibacilus borealis"
for record in SeqIO.parse(handle, "fasta"):
if organism_name in record.id:
# Do whatever you want with these
print record.id
print record.seq
handle.close()
```

There's more info here: http://biopython.org/wiki/SeqIO
Question: Hi everyone,

I was wondering if anyone if familiar of any annotation term in the human genes annotation e.g. from gencode or ensembl to be able to extract **pseudogenes** and separate them from non-pseudogenes.

Any thoughts?

Thanks in advance,
Sergio

Answer: You can use Ensembl BioMart with the following query:

[very long BioMart query..... Modify parameters as you like][1]

There are many sub-types of pseudo genes, the query outputs the gene type in the last column.

[1]: http://www.ensembl.org/biomart/martview/d5e8623bf33a920fd3fdf6149c374e2a?VIRTUALSCHEMANAME=default&ATTRIBUTES=hsapiens_gene_ensembl.default.feature_page.ensembl_gene_id%7Chsapiens_gene_ensembl.default.feature_page.ensembl_gene_id_version%7Chsapiens_gene_ensembl.default.feature_page.ensembl_transcript_id%7Chsapiens_gene_ensembl.default.feature_page.ensembl_transcript_id_version%7Chsapiens_gene_ensembl.default.feature_page.start_position%7Chsapiens_gene_ensembl.default.feature_page.end_position%7Chsapiens_gene_ensembl.default.feature_page.strand%7Chsapiens_gene_ensembl.default.feature_page.chromosome_name%7Chsapiens_gene_ensembl.default.feature_page.description%7Chsapiens_gene_ensembl.default.feature_page.gene_biotype&FILTERS=hsapiens_gene_ensembl.default.filters.biotype.%22IG_C_pseudogene,IG_J_pseudogene,IG_V_pseudogene,polymorphic_pseudogene,processed_pseudogene,rRNA_pseudogene,transcribed_processed_pseudogene,transcribed_unitary_pseudogene,transcribed_unprocessed_pseudogene,translated_processed_pseudogene,translated_unprocessed_pseudogene,TR_J_pseudogene,TR_V_pseudogene,unitary_pseudogene,unprocessed_pseudogene%22&VISIBLEPANEL=resultspanel
Question: I came across this website infinome.me which relates SNPs to a few different diseases. On their site, if you click on a pubmed article, they give you a list of all the SNPs mentioned in that article along with the p-value (relating the SNP to the disease). Here's a screenshot:

^^ scroll to the right to see full image.

Is there an online database where I can download this kind of information, like for example if I'm interested in schizophrenia, I can get some articles with the list of SNPs they tested and their corresponding p-values?

Answer: Try [Grasp](http://grasp.nhlbi.nih.gov/Search.aspx) or [GWAS central](http://www.gwascentral.org/).
Question: Hello!

I apologize in advance if this is not something that should go to the knowledge base, but I keep hitting a similar wall in analyzing FASTQ files on SRA that are open to public. When dealing with exome data I keep having trouble finding the information which capture kit is used, especially if the paper is hidden behind a subscription (Which it often is). When I do find out which capture kit was used, I have trouble finding the BED file.

For example, I found out a project used Qiagen Qiaamp DNA Mini kit. I registered with the website and went all over the internet, but I haven't found any BED files, or any resources that explain their acquisition. Customer support is ignoring me, possibly because it is a silly question. In general, how does one find these BED files and why are they not public?

Thanks in advance!

Answer: The Qiaamp is a DNA purification kit, it does not do any target enrichment. Therefore, customer service is probably ignoring you, as the question does not make sense. I would check the GEO entry of the respective SRA if there are any information given in the methods section. Still, the methods section in the publication must provide this information. Why are many journals not public...well because capitalism. How many publications would you need to check?

EDIT: As pointed out below, I was mixing up GEO and SRA. It is not GEO to check for details, but rather the [SRA portal.][1] GEO does [not accept][2] WES/WGS data, but some DNA samples such as ChIP-seq.

[1]: https://www.ncbi.nlm.nih.gov/sra
[2]: https://www.ncbi.nlm.nih.gov/geo/info/seq.html#categories
Question: I am developing an R package and I plan to eventually upload it to BioConductor. However, after I went over the guidelines, I realized that it seems to be more restricted comparing with the requirement from CRAN. I followed the book: "[R Package][1]" by Hadley Wickham, but I believe this book is designed to help people submit CRAN packages, not too specifically for BioConductor.

Therefore I wonder if there is any special issue I should pay additional attention here. Especially formats I guess.

Thank you.

[1]: http://r-pkgs.had.co.nz/

Answer: I have a package at Bioconductor and I really appreciated the support during the initial review, where I was asked some improvements of quality that were not too demanding, but still made me learn more about S4 classes. I find the 6-month release schedule to be a good regular reminder, and in the worst case there is now [a procedure to retire packages](https://www.bioconductor.org/developers/package-end-of-life/). New packages can be [submitted through GitHub](https://github.com/Bioconductor/Contributions/blob/master/CONTRIBUTING.md), and it is possible to [push updates via a nice Git-Svn bridge](https://www.bioconductor.org/developers/how-to/git-mirrors/). Lastly, the use of Bioconductor's base S4 classes is not mandatory (but of course comes naturally once one gets more experienced). Altogether, Bioconductor is a brilliant platform and it is worth the effort aiming for it.
Question: I am searching for full mitochondrial genome of a species of [Malvaceae][1] family. But I can not locate where it is in gene bank. Searching in google, I found [this link][2], but here I did not found anything for plant. Can you help me?

To be more specific, I need the following scientific classification holder:

```
Kingdom: Plantae
(unranked): Angiosperms
(unranked): Eudicots
(unranked): Rosids
Order: Malvales
Family: Malvaceae
Subfamily: Grewioideae
Genus: Corchorus L.
```

from [Wikipedia][3]

[1]: https://en.wikipedia.org/wiki/Malvaceae
[2]: http://www.ncbi.nlm.nih.gov/genomes/ORGANELLES/mztax_short.html
[3]: https://en.wikipedia.org/wiki/Corchorus

Answer: Mitochondrial (and plastid) genomes are hosted at [NCBI Organelle Genome Resources][1]. Just browse Malvaceae and filter for mitochondria. I came up with two complete genomes.

Alternatively, you can use these queries:

- http://www.ncbi.nlm.nih.gov/genomes/GenomesGroup.cgi?opt=organelle&taxid=3629 (view)
- http://www.ncbi.nlm.nih.gov/genomes/GenomesGroup.cgi?opt=organelle&taxid=3633&cmd=download1 (directly download accessions)

[1]: http://www.ncbi.nlm.nih.gov/genome/organelle/
Question: Hello,

I have 600+ bacterial isolates (*fastq.gz files) that I want to assemble.

I have a script which uses a for loop to trim adapters and run SPAdes on each isolate, but this takes a very long time, and my computer should be powerful enough to run the script for multiple isolates at one time, which is what I want to do.

I have looked at GNU parallel which seems to work well and be faster than the for loop. As an example, I do the trimming like this:

parallel 'trimmomatic PE {}R1*.f*q.gz {}R2*.f*q.gz {}pair_R1.fq.gz {}unpair_R1.fq.gz {}pair_R2.fq.gz {}unpair_R2.fq.gz ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36' ::: $(ls *.fastq.gz | rev | cut -c 16- | rev | uniq)

My question is; When running parallel in a folder containing my 600+ filepairs, will parallel try to run the program on all 600+ isolates simultaneously or will it limit the amount of files run at one time based on how much the computer can manage at one time? Or is there a way to limit how many files parallel should work on at one time, apart from specifying the specific files?

Thank you!

Answer: Parallel will run as many as it can concurrently* (i.e the number of cores you have, unless you specify otherwise with `-j`. There is also a "max load" flag you can use in the manual to manage the load on the system.

*Edit, the caveat of this being that you're running/invoking enough commands. Obviously if you have 32 cores, but only 16 files to work on, only 16 threads will be spawned.
Question: Hello,

I would like to check whether L1 repetitive elements are modulated between my treatment and control via RNA-Seq. I have read several papers that have done so but their methods are not clear enough for me as a biologist to reproduce. I have analyzed my data using the Tuxedo suite and have analyzed the "unique" genes. I am wondering what modifications have to be taken into account to accommodate the repetitive nature of LINEs.
1- I have an understanding that some aligners filter out reads that map to several places in the genome. Are my LINE reads being filtered out by tophat?
2- If so, how do I align them?
3- when using cufflinks, intead of using RefSeq, I am assuming I would have to use a repetitive element model?

Thank you!

Answer: You might get some ideas from a solution described in a paper from Peter Park's lab, Estimating enrichment of repetitive elements from high-throughput sequence data which has an online tool available, with source code (Repeat Enrichment Estimator). It appears to be for ChIP-seq though; not sure how adaptable it would be for RNA-seq.

Edit on Apr 27 2015:

I recently had to revisit this problem and found a useful tool that didn't exist at the time of the original answer:

RepEnrich (paper, github)
Question: I have 9 samples for bulk RNA-seq analysis with DESeq2. I have 2 conditions and 2 batches. Below is the sample table.

![enter image description here][1]

My design is based on both batch and condition and the PCA plot shows that C1,C2, and C3 samples are separated. below is the code and the PCA plot.

dds <- DESeqDataSetFromTximport(txi.rsem, sampleTable, ~Batch+condition)
dds <- DESeq(dds)
vsd <- vst(dds, blind=FALSE)
plotPCA(vsd, intgroup=c("condition", "Batch"))+ geom_text(aes(label=name),vjust=2)

![enter image description here][2]

I know that it is recommended not to remove the batch effect but I was trying to check and see what would be the results after removing the batch effect. I again see that the control samples C1,C2,C3 are not clustered with the other control C4 from another batch. Below is the code and PCA plot after removing the batch effect.

mat <- assay(vsd)
mat2 <- limma::removeBatchEffect(mat, vsd$Batch)
assay(vsd) <- mat2
plotPCA(vsd, intgroup=c("condition", "Batch")) + geom_text(aes(label=name),vjust=2)

![enter image description here][3]

My question is does it make sense to go further and do the DEG analysis? or my samples are not clustered/separated rationally and the results would not make any sense?
I think I cannot remove C4 because that is the only connection between the samples of different batches (if that makes sense!). I was hoping that C1,C2, and C3 would cluster together with C4 and then the comparison between the Control and Treatment of different batches would be accurate.

[1]: https://i.imgur.com/J8c9BY9.png
[2]: https://i.imgur.com/iBS6kg0.png
[3]: https://i.imgur.com/ice9NxP.png

Answer: Batch correction requires replicates of each group in both batches. Your experiment is fully confounded with all treatments in B2. Sorry to say but this is simply an unfortunate experimental design and I do not think any in silico methods can save you from that.