cannin/biostars_qa huggingface dataset:
Question: Dear All,

I am using BAM files for chip-seq analysis. The chromosome notation in a usual BAM file is like: chr1. In my file the chromosomme notation is 1. Is there a way to change that in the BAM file? For further analysis it is very important to change the notation.

Many thanks!

Greetz Lisanne

Answer: The following line works well for me. It will replace the names of chromosome 1 to 22, X, Y and MT (mitochondrial DNA). It will create a new file with the name `test_chr.bam` instead of `test.bam`.

use Samtools

samtools view -H test.bam |\
sed -e 's/SN:1/SN:chr1/' | sed -e 's/SN:2/SN:chr2/' | \
sed -e 's/SN:3/SN:chr3/' | sed -e 's/SN:4/SN:chr4/' | \
sed -e 's/SN:5/SN:chr5/' | sed -e 's/SN:6/SN:chr6/' | \
sed -e 's/SN:7/SN:chr7/' | sed -e 's/SN:8/SN:chr8/' | \
sed -e 's/SN:9/SN:chr9/' | sed -e 's/SN:10/SN:chr10/' | \
sed -e 's/SN:11/SN:chr11/' | sed -e 's/SN:12/SN:chr12/' | \
sed -e 's/SN:13/SN:chr13/' | sed -e 's/SN:14/SN:chr14/' | \
sed -e 's/SN:15/SN:chr15/' | sed -e 's/SN:16/SN:chr16/' | \
sed -e 's/SN:17/SN:chr17/' | sed -e 's/SN:18/SN:chr18/' | \
sed -e 's/SN:19/SN:chr19/' | sed -e 's/SN:20/SN:chr20/' | \
sed -e 's/SN:21/SN:chr21/' | sed -e 's/SN:22/SN:chr22/' | \
sed -e 's/SN:X/SN:chrX/' | sed -e 's/SN:Y/SN:chrY/' | \
sed -e 's/SN:MT/SN:chrM/' | samtools reheader - test.bam > test_chr.bam

If you want to do this for a bunch of bam files at once, use a loop. The variable `$filename` is your file without the extension (assuming there is no period except the one in ".bam").

for file in *.bam
do
filename=`echo $file | cut -d "." -f 1`
samtools view -H $file | \
sed -e 's/SN:1/SN:chr1/' | sed -e 's/SN:2/SN:chr2/' | \
sed -e 's/SN:3/SN:chr3/' | sed -e 's/SN:4/SN:chr4/' | \
sed -e 's/SN:5/SN:chr5/' | sed -e 's/SN:6/SN:chr6/' | \
sed -e 's/SN:7/SN:chr7/' | sed -e 's/SN:8/SN:chr8/' | \
sed -e 's/SN:9/SN:chr9/' | sed -e 's/SN:10/SN:chr10/' | \
sed -e 's/SN:11/SN:chr11/' | sed -e 's/SN:12/SN:chr12/' | \
sed -e 's/SN:13/SN:chr13/' | sed -e 's/SN:14/SN:chr14/' | \
sed -e 's/SN:15/SN:chr15/' | sed -e 's/SN:16/SN:chr16/' | \
sed -e 's/SN:17/SN:chr17/' | sed -e 's/SN:18/SN:chr18/' | \
sed -e 's/SN:19/SN:chr19/' | sed -e 's/SN:20/SN:chr20/' | \
sed -e 's/SN:21/SN:chr21/' | sed -e 's/SN:22/SN:chr22/' | \
sed -e 's/SN:X/SN:chrX/' | sed -e 's/SN:Y/SN:chrY/' | \
sed -e 's/SN:MT/SN:chrM/' | samtools reheader - $file > ${filename}_chr.bam
done
Question: Hi, I used samtools to generate a mpileup file. In the base column, I notice some "!". I sought out the sigification in domcumentation/already asked question and haven't find their meaning. Does anybody know?

Exemple:

FNXXX	2595	T	1	^!a	J

Answer: Try to add ```-Bx``` to your command.

Short story long:

mpileup is originally made to call variants. To not count the overlapping bases between the forward and reverse reads twice samtools changes the quality score to "!". This "!" has a phred score of 0 and will be ignored in the variant calling process.

Here more info:
https://sourceforge.net/p/samtools/mailman/message/32793789/
Question: Hello all,

I was a little confused by this plink output im producing. It was my understanding that --indep-pairwise analysis using a cutoff value for r2 would identify all the SNPs in a .ped file with an r2 correlation value of >.7

$ ./plink --file input.plink --indep-pairwise 225 23 .7 --out plinkoutput1

SNPs with a r2 >.7 would have their IDs saved to a file 'prune.in' while those with a r2 < .7 would be saved to a second file called 'prune.out'. Now, if you want to get the tagging SNPs and a list of all the SNPs they tag, you run something along the lines of

$ plink --file input.plink --show-tags my_tags.txt --list-all --tag-r2 .7 --out mytags_and_what_they_tag

So it makes sense that you could use the prune.in list of SNPs for the list of tags right?

$ plink --file input.plink --show-tags input.prune.in --list-all --tag-r2 .7 --out mytags_and_what_they_tag

and in the output every SNP should appear as either a tag for something, or something that is being tagged by another snp. But when I look at the file:

SNP CHR BP NTAG LEFT RIGHT KBSPAN TAGS
rs79687004 2 5501148 0 5501148 5501148 0 NONE
rs1404257 2 5503200 1 5498227 5503200 4.973 rs10929518
rs61268182 2 5504924 2 5491746 5504924 13.178 rs72776977|rs955146
rs17356301 2 5506274 0 5506274 5506274 0 NONE

I see things like rs79687004. In its 'TAGS' column it has NONE. So its not serving as a tagging SNP for anything. Fair enough, but when I search for it in the file this is the ONLY instance of 'rs79687004'. Its not being listed in some other SNP's TAGS column either. Heres where I am confused: if everything in prune.in is a SNP that correlated to something else at an r2 > .7, wouldnt every single SNP either be identified as tagging something or tagged by something in this output? Ie shouldnt a SNP that is not tagging anything, in addition to having its on line in the output ALSO appear in some other line as a tagged SNP?

The only thing I can think as an explanation is that my scanning window used during --pairwise-indep was 225 SNPs, and this file was a collection of some distant loci in a single chromosome, so maybe there were instances where the --indep-pairwise window fell over a VERY distantly located set of SNPs and found they were correlated, so they went into the 'prune.in' file. But then, during the --show-tags analysis there is a 240kb default 'window' where SNPs further apart than that are not examined for tagging. but I find it hard to swallow that I would get so many lines where a SNP is not tagging OR tagged just because some comparisons won't be made the second time around in --show-tags.

What am I missing here? I've seen other peoples output from the same process posted with "None" in the TAGS column, but no one seems to bat an eye. I feel like I am overlooking something obvious here and I was hoping you could point it out to me.

Many thanks in advance!

Answer: It's the .prune.out SNPs which are guaranteed to have r^2 > .7 with another SNP here. A .prune.in SNPs might tag other SNPs, or it might just tag itself.
Question: Hi,

I have data with 4 levels group: condition1, treat1, condition2, treat2.

I am trying to use DESeq2 to identify differentially expressed genes between two levels: treat1 vs condition1.

I tried two methods. First, I only use subset of all data, condition1 and treat1 to get the differentially expressed genes between treat1 vs condition1. Second, I use all the data including condition1, treat1, condition2 and treat2; but I extract the differentially expressed genes between treat1 vs condition1 by `results(dds, contrast("group","treat1","condition1"))`.

Interestingly, the identified differentially expressed genes are completely different. When I checked size factors for each sample in these two different methods, they are also completely different. Is this an expected behavior for DESeq2 pipeline? If so, which method should I use?

Thank you.

Answer: This is not surprising, because the model estimates (size factors, dispersion, priors, etc) use all samples. I suggest exploring a PCA plot to see the relationships between samples (see vignette).
Question: I've been experimenting with the [GenomicDataCommons][1] package to handle query work against the GDC API. For some reason there is an issue with timeouts for requests of a certain length. There doesn't seem to be a direct way around this using the piped syntax. Anyone else have luck with this?

There is a results(size = n) method, its syntax seems to allow only the first n records to be accessed.

Here is an example query (should return 400-500 records):

proj <- 'TCGA-COAD'
case_data <- cases() %>%
GenomicDataCommons::filter(~ project.project_id == proj) %>%
GenomicDataCommons::expand('diagnoses') %>%
results_all() %>%
as_tibble()

Gives, after a few moments:

Error in is.response(x) : Internal Server Error (HTTP 500).

[1]: https://bioconductor.org/packages/release/bioc/html/GenomicDataCommons.html

Answer: I should clean up the documentation, but `results_all()` is a convenience wrapper that is not too smart in that it simply tries to return all results in one trip to the server. This can fail for multiple reasons related to the size of result sets. The better approach (and the only one in the case of large results sets) is to page through the results:

proj <- 'TCGA-COAD'
query = cases() %>%
GenomicDataCommons::filter(~ project.project_id == proj) %>%
GenomicDataCommons::expand('diagnoses')
count = query %>% count()
size = 50
reslist = lapply(seq(1,count, size), function(start) {
query %>%
results(size=size, from = start) %>%
as_tibble()
})
case_data = bind_rows(reslist)

Unfortunately, the `size` parameter really requires trial-and-error to find the largest "working" setting since the results can vary quite significantly in volume. Instead, I usually just choose a smallish number like 50 or so and wait a few extra seconds. These calls can, in theory, be parallelized using something like BiocParallel to get really fancy (and introduce complexity).
Question: Hi All,

After installing kraken am trying to build the database as specified in the manaul but getting the following messages.
Any inputs on this??

/Tools/kraken-master/KRAKEN$ ./kraken-build --standard --threads 40 --db /home/karthic/Databases/KRAKEN
Found jellyfish v1.1.11
Step 1/3: performing rsync dry run...
Rsync dry run complete, removing any non-existent files from manifest.
Step 2/3: Performing rsync file transfer of requested files
rsync: failed to connect to ftp.ncbi.nlm.nih.gov (165.112.9.229): Connection timed out (110)
rsync: failed to connect to ftp.ncbi.nlm.nih.gov (2607:f220:41e:250::7): Network is unreachable (101)
rsync error: error in socket IO (code 10) at clientserver.c(128) [Receiver=3.1.1]
rsync_from_ncbi.pl: rsync error, exited with code 10

Thanks in Advance,
KK

Answer: Since NCBI updated their FTP website and decided to phase-out Genbank Identifiers (GIs), the default Kraken database update scripts do not work.

My colleague @Sej Modha has written a python script that helps with updating the kraken databases:
http://bioinformatics.cvr.ac.uk/blog/update-kraken-databases/
Question: where can I get the parameters used by UCSC genomes to run RepeatMasker for dm6? I searched FAQ but didn't find the answer.

Answer: [Here][1].

[1]: https://genome.ucsc.edu/FAQ/FAQdownloads.html#download16
Question: Im using the heatmap.2 function in R to create a clustered heatmap of some microarray data. My commands are written as this:

heatmap.2 (as.matrix(highly_variable_genes,col=rev(morecols(50)),trace="none", main="Top 500 most variable genes across samples",ColSideColors=col.cell,scale="row"))

But even though I have specified that I don't want that hideous cyan squiggly 'trace' over my cells, which obscures everything, it keeps coming up. I've tried re-writing the command a few different ways and cannot get rid of it. Does anyone see what I have done wrong?

This command is taken wholesale from a tutorial that has no trace over its heatmap cells so I really dont know whats going on here...unless the commands have been updated in R.

Many thanks for your help.

Answer: I agree that the tracing leaves much to be desired and I have never used it.

There is something peculiar about the way you have written your code - take a closer look at it. You are passing all of the following to
`as.matrix()`, which fails to issue any warning: `x,col=rev(morecols(50)),trace="none", main="Top 500 most variable genes across samples",,scale="row"`

Here is a simple solution:

require(gplots)

# generate random data

x <- matrix(rexp(200, rate=.1), ncol=20)

# reproduce the 'bug' with trace lines

heatmap.2(as.matrix(x,col=rev(morecols(50)), trace="none", main="Top 500 most variable genes across samples",,scale="row"))

# get rid of trace lines

heatmap.2(as.matrix(x), trace="none", main="Top 500 most variable genes across samples", scale="row")
Question: Hi,

I am currently curating a list of prognostic genes and gene variants responsible for prostate and/or breast cancer. I am currently concentrating only upon germline variants, but am also collecting the somatic variants as when I run into them. Is there a database which already has such a list of genes and also the papers relating to them? The information which I am looking for and curating are the genes, gene locations, mutation type and the link to the paper. Is there any publicly available database which lists all this information.

Thanks,
Venkat

Answer: You can look at some of the resources from [Open Targets][1], such as Open Targets Genetics and Open Targets Platform. With the latter, you can find genes associated with [prostate carcinoma][2] and [breast carcinoma][3] based on [somatic mutations][4], germline variants (giving rise to what we call [genetic associations][5]) and many other types of different [data sources][6]. The main entry points here is a gene or disease name.

With [Open Targets Genetics][7], you can find common (germline) variants linked to those traits of interest and the genes that are (potentially) functionally implicated by these common variants. Have a look at these two random examples of studies: 1) [Breast cancer][8] by Michailidou K (2017) in Nature and 2) [Prostate cancer][9] by Schumacher FR (2018) in Nat Genet. And the cool thing is that you can also compare different studies, such as the two above in the "[Compare overlapping studies][10]" page. The main entry points here is the variant (rs ID or GRCh38 coordinates), gene or study/trait names.

Both Open Targets Platform and Open Targets Genetics complement each other (with links from the latter to the former) and there is plenty of things you can do with them, some of them depicted in our short tutorial videos on [YouTube][11].

[1]: https://www.biostars.org/p/196813/
[2]: https://www.targetvalidation.org/disease/EFO_0001663/associations
[3]: https://www.targetvalidation.org/disease/EFO_0000305/associations
[4]: https://docs.targetvalidation.org/data-sources/somatic-mutations
[5]: https://docs.targetvalidation.org/data-sources/genetic-associations
[6]: https://docs.targetvalidation.org/data-sources/data-sources
[7]: https://www.biostars.org/p/350961/
[8]: https://genetics.opentargets.org/study/GCST004988
[9]: https://genetics.opentargets.org/study/GCST006085
[10]: http://%20https://genetics.opentargets.org/study-comparison/GCST004988?studyIds=GCST006085
[11]: https://www.youtube.com/channel/UCLMrondxbT0DIGx5nGOSYOQ
Question: Dear all,

This is driving me crazy, maybe someone has an idea.

I simply want to print the rows of my dataframe with

write_xlsx(dd, "Z:/dd.xlsx",row.names=TRUE)

But this error keeps coming:

"unused argument (row.names = TRUE)"

Why is it not working?

It works if I remove "rownames=TRUE", but then it doesnt print the rownames.

Thank you so much!
Bine

Answer: It appears write_xlsx doesn't have a parameter `row.names`. Why don't you just add a column for row names to the dataframe `dd` before writing it? Like:

dd$row_names <- rownames(dd)
dd <- dd[, c(ncol(dd), 1:(ncol(dd)-1))] # Make row_names first column
Question: Hi,
I have a large paired-end dataset in the BAM format and a list of read IDs which belong to a single mate of a pair. What I want to do is to extract their second mates from the whole dataset. Could you please advise me some efficient ways to do this like using, let's say, Bio-SamTools or something like that? Something memory- and time-efficient.
Thanks!

Answer: **2018**: use picard https://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

____

The following C++ code should filters only print the reads contained in your file (I've added to my variation toolkit http://code.google.com/p/variationtoolkit )

http://code.google.com/p/variationtoolkit/source/browse/trunk/src/bamgrepreads.cpp

Compilation:

g++ -O3 -Wall -I path/to/samtool -L path/to/samtool bamgrepreads.cpp -lbam -lz

Usage:

./a.out -R file_containing_the_reads_name.txt (stdin|bam1 bam2...)

Options:

-f INT required flag
-F INT filtering flag
-R FILE reads file
-e only one match per name (goes faster)
Question: I have a data.frame where I need to collapse rows by sample names in the indiv.ID column. Here is an excerpt:

indiv.ID 86912632 86920881 86922082 86927699
1 Alxis_3702 CTGA
2 Alxis_3702 TCTG
3 Alxis_3702 G
4 Alxis_3702 C
5 Alxis_3702
6 Alxis_3702
7 Alxis_3702
8 Alxis_3702
9 Alxis_3702
10 Alxis_3702

In theory it should be easy to merge all rows with the same indiv.ID, but in the case of rows 1 and 2, I also need a new column (because this locus is heterozygous), so that each row only has one individual, like so:

indiv.ID 86912632 86912632.1 86920881 86922082 86927699
Alxis_3702 CTGA TCTG G C

What is the best way to do this?

Answer: We can `gather` then `spread`, see example (I changed the input data to demonstrate that this can be used for data with multiple *indiv.IDs*):

# example data, with 2 IDs
df1 <- read.table(text = " indiv.ID 86912632 86920881 86922082 86927699
1 Alxis_3702 CTGA NA NA NA
2 Alxis_3702 TCTG NA NA NA
3 Alxis_3702 NA G NA NA
4 Alxis_3702 NA NA C NA
5 Alxis_3702 NA NA NA NA
6 Alxis_3702 NA NA NA NA
7 Alxis_3703 AA NA NA NA
8 Alxis_3703 NA TT CC NA
9 Alxis_3703 NA NA NA NA
10 Alxis_3703 NA NA NA GG")

library(dplyr)
library(tidyr)

gather(df1, key = "k", value = "v", -indiv.ID) %>%
filter(!is.na(v)) %>%
group_by( indiv.ID) %>%
mutate(k = make.unique(k)) %>%
spread(key = k, value = v) %>%
ungroup()

# # A tibble: 2 x 6
# indiv.ID X86912632 X86912632.1 X86920881 X86922082 X86927699
#
# 1 Alxis_3702 CTGA TCTG G C NA
# 2 Alxis_3703 AA NA TT CC GG
Question: I'm trying to infer the the type of experiment for my paired-end RNAseq data.
I have employed STAR to index the genome and map the paired reads (as unstranded) to derive the experiment. Depending on the bed file I employ to infer the experiment I get different results.
As genome I'm employing the genecode:

> GRCm39.primary_assembly.genome.fa

Anyone can guess why are these differences happening?

**Method 1 (bed converted from genecode gtf)**

Convert gtf to bed

bedparse gtf2bed gencode.vM27.annotation.gtf > gencode.vM27.annotation.gtf.bed

**Code**

infer_experiment.py -r /local/ref/gencode.vM27.annotation.gtf.bed -i /local/out/star/ADU2124Aligned.sortedByCoord.out.bam

**Result:**

This is PairEnd Data
Fraction of reads failed to determine: 0.1017
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0080
Fraction of reads explained by "1+-,1-+,2++,2--": 0.8903

**Method (bed file provided by rseqc)**

**Code**

infer_experiment.py -r /local/ref/mm10_Gencode_VM18.bed -i /local/out/star/ADU2124Aligned.sortedByCoord.out.bam

**Result**:

This is PairEnd Data
Fraction of reads failed to determine: 0.0577
Fraction of reads explained by "1++,1--,2+-,2-+": 0.3863
Fraction of reads explained by "1+-,1-+,2++,2--": 0.5560

Thanks!

Answer: Well, it seems that I was employing the wrong bed from rseqc and I didn't notice until now...

> **BED:** mm10_Gencode_VM18.bed
>
> **Reference:** GRCm39.primary_assembly.genome.fa
Question: Here is an example of volcano plot produced using edgeR results on DEG analysis. So my question is why there is no values in the middle of plot (indicated by red circle)?

![Volcano plot][1]

[1]: https://i.pinimg.com/originals/07/69/3d/07693dfcb3c8fe55fcae39d7caa64239.jpg

Answer: You are plotting the fold change of differential expression versus the confidence in detecting this differential expression. (q-value) Genes which have a low fold change will generally have a higher q-score.
Question: I want to know if there is a way to keep sites with 0% missing a VCF file.

On vcftools I can do:

`vcftools --vcf input_file.vcf --max-missing 1.0 --out output_noMissing`

Is there a way to do this with BCFTOOLS as it is much faster?

Answer: $ bcftools view -e 'GT[*]="mis"' input.vcf

This will exclude (`-e`) sites where any genotype (`GT[*]`) is missing (`="mis"`).
Question: Has anybody made a systematic survey on quality control parameters in published NGS-data?
I am especially interested in knowing about duplication levels people see in ChIP-Seq and RNA-Seq data.

Is there anything like this available on the encode data, or does one have to download that and QC oneself?

Answer: This would be a very interesting, but challenging study. It seems that in the early days, every NGS paper that was published would receive a lot of probing from reviewers about library depth, coverage, percent duplicates, library complexity, quality distribution, etc. But now it is common to exclude mention of these kinds of metrics as if the methods were well-developed and standardized. But, in my opinion this is far from the case. Its not like with an Affy chip where you have a well-defined system, with well-developed processing methods and QC protocols. I think that RNA-seq is fundamentally better than arrays as a data type but the analysis and QC of such data still resides very much in the wild-west. One group to watch in this area is MAQC-III (also known as SEQC). They are "assessing the technical performance of next-generation sequencing platforms by generating benchmark datasets with reference samples and evaluating advantages and limitations of various bioinformatics strategies in RNA and DNA analyses". They are looking at multiple platforms, performing spike-in experiments, comparing performance to benchmarks (e.g., arrays), etc and have about 25 analysis teams involved. I think we can probably expect a special issue on the results of their first round some time in late 2012.
Question: Is anyone aware of a publication or pre-publication data set that includes the following for a single human tumor/normal sample pair:

Illumina whole genome sequence (WGS) to at least 30-50X depth for both tumor and normal (e.g. blood). Preferably this would be HiSeq paired end data generated with v3 chemistry

Illumina exome sequence data for the same tumor/normal pair but to considerably higher depth (say at least 150-200X)

Illumina RNA-seq data generated for the same tumor sample. Bonus points if there is also RNA-seq for a matched normal sample (same tissue as the tumor as opposed to blood DNA that would typically be used for a normal comparison at the DNA level).

Since this data would primarily be used for methods development, the tumor/normal pair could be from a tumor cell line and matched lymphoblastoid 'normal' cell line derived from the same individual. Some data along these lines is being made available here:

TCGA Mutation/Variation Calling Benchmark 4 at CGHub

https://cghub.ucsc.edu/benchmark_download.html

However, this is whole genome data only. No exome or RNA-seq data yet. Plenty of RNA-seq data can be found in GEO and elsewhere but I'm not aware of any projects where there is corresponding WGS or Exome data. Plenty of exome data and WGS data are being generated for TCGA but again I'm not aware of any publications describing combinations of all three types.

Large scale cancer sequencing projects that might have performed such a comparison:

The Cancer Genome Atlas (TCGA)

Cancer Genome Project (CGP)

International Cancer Genome Consortium (ICGC)

Answer: We eventually created such a data set ourselves and made it publicly available via FTP [here](https://github.com/genome/gms/wiki/HCC1395-WGS-Exome-RNA-Seq-Data)

The data are from a matched tumor/'normal' pair of cell lines: HCC1395 and HCC1395/BL whole genome (WGS), exome, and/or RNA-seq data. All data are 2x100 bp reads generated on an Illumina HiSeq 2000 instrument. The exome data was generated by use of a NimbleGen SeqCap EZ Human Exome Library v3.0 reagent.

If you find this data useful, please cite:

Griffith et al. Genome Modeling System: A Knowledge Management Platform for Genomics. PLoS Comput Biol. 2015 Jul 9;11(7):e1004274. doi: 10.1371/journal.pcbi.1004274. eCollection 2015 Jul. PubMed PMID: 26158448;

Since these data corresponds to cell line material we were able to make them available without using dbGaP.