cannin/biostars_qa huggingface dataset:
Question: I have a table like below, it is a bed file of genome coordinate, I would like to keep only rows with numbers.

Input:

1	141009669	141009952
9	141016322	141016973
GL000195.1	81719	82468
GL000195.1	142613	142923
GL000220.1	119445	119746
HG115_PATCH 101957832	101958132
HG1308_PATCH 130205069	130205369
HG1308_PATCH 130205406	130205773
HG748_PATCH 77577953	77578264
X	200983	202660
y 205180	205702

output:

1	141009669	141009952
9	141016322	141016973
X	200983	202660
y	205180	205702

Thanks in advance!

Answer: A straight forward way:

mainChr = c(as.character(1:22),'x','X','y','Y')
data = read.delim('your.bed',stringsAsFactor = F,header = F)
data_with_mainChr = data[data$V1 %in% mainChr,]

if you want to use `readr` and `dplyr` which is more efficient when dealing with big files:

mainChr = c(as.character(1:22),'x','X','y','Y')
library(dplyr); library(readr)
data = read_tsv('your_bed_file',col_names = F)
data_with_mainChr = dplyr::filter(data, X1 %in% mainChr)
Question: Hi all,

I have [Arlequin][1] output files (*.arp) from [fastSimcoal2][2] that I'm trying to convert to genpop files (to read into [adegenet][3]). The Arlequin files are using the "DNA" marker (50 bp long), and each includes ~2,500 samples.

When I try to use the program [PGDspider][4] to convert these files (either the GUI or the CLI, v2.1.1.5), things proceed smoothly until the 1,000th allele is reached, where the program generates an error:

```
INFO 07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596) -Allele "TGGCATCAGCTTCTTCTAAGCAGGCCCGCTTTGTCTAGCGACAGATCTCG" converted to "999"
INFO 07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596) -Allele "ATCTTGTGGTCATGGAGCACCCAGAGCCAAACTCCTAAAGAAGACCCCGC" converted to "1000"
ERROR 07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:611) -data "1000" in ind "2_375" (population: Sample 2) is too long!!!
INFO 07-03 11:18:10,631 (GenepopWriter.java:parseAlignedSampleData:618) -In GENEPOP, alleles cannot be coded with more than 3 digits.
```
Is there a way around this? I've processed RADseq data with many 1,000s of alleles using adegenet before, so I know these types of files can be written to genpop--I just don't know how to specify this to PGDspider. There's a field in the .spid file that I might be able to use to fix this (see below), but I don't know what to populate that field with.
```
# Specify the locus/locus combination you want to write to the GENEPOP file:
GENEPOP_WRITER_LOCUS_COMBINATION_QUESTION=
```

Alternatively, if there are other suggestions for performing this conversion without using PGDspider, I'm all ears.

Thanks for your help!

[1]: http://cmpg.unibe.ch/software/arlequin35/
[2]: http://cmpg.unibe.ch/software/fastsimcoal27/
[3]: https://cran.r-project.org/web/packages/adegenet/index.html
[4]: http://www.cmpg.unibe.ch/software/PGDSpider/#Download_and_Installation_Instructions

Answer: Answer: I've been able to do this using the R package [strataG][1] (no longer on CRAN, but can be obtained on GitHub). Specifically, the steps are:

1. Read in .arp file, e.g. arpObject <- arlequinRead(file.arp)
2. Convert arlequin object to gtype object (arp2gtypes function)
3. Convert gtype object to genind object (gtypes2genind function)

[1]: https://github.com/EricArcher/strataG
Question: I have some profiles with scores in BED format and I would like to convert them into bigWig format to visualise these data in jBrowse.

```
9 35084375 35084376 ENSDART00000000004 160.9480969 -
9 35060460 35060461 ENSDART00000000005 752.3776435 +
4 14148184 14148185 ENSDART00000000019 5226.62898 -
2 50770137 50770138 ENSDART00000000042 3513.471655 -
```

Typically (genomeCoverageBed + bedGraphToBigWig), the intervals from BED file are summed up into bigWig, but this is not what I want. I would like the scores from BED to be visible directly in bigWig track.

Can you recommend me some way of converting such BED in bigWig?

Answer: You can use awk to get a bedgraph:

awk '{printf "%s\t%d\t%d\t%2.3f\n" , $1,$2,$3,$5}' myBed.bed > myFile.bedgraph

In case the BED file was not sorted properly:

sort -k1,1 -k2,2n myFile.bedgraph > myFile_sorted.bedgraph

Finally, use the UCSC bedGraphToBigWig tool:

bedGraphToBigWig myFile_sorted.bedgraph myChrom.sizes myBigWig.bw

Beware of the naming conventions (sometimes you need to add a chr to your chromosomes' name).

Cheers,
Michael

[edit] updated awk-command after comment
Question: Hi all,

I have a quick question.

When I run bwa mem for oxford nanopore with the option:

bwa mem -x ont2d

Should I include the option:

bwa mem -x ont2d -M

For split reads? If I want to run Lumpy on the alignments would I need to supply "-M"?

Answer: According to the [examples of Lumpy][1] you don't have to use `-M`:

bwa mem -R "@RG\tID:id\tSM:sample\tLB:lib" human_g1k_v37.fasta sample.1.fq sample.2.fq \
| samtools view -S -b - \
> sample.bam

[1]: https://github.com/arq5x/lumpy-sv
Question: Hello, I would like to know how to delete my biostars account. Thank you in advance.

Answer: Write an email to admin@biostars.org with the request and your biostars account details
Question: Other than Limma being designed for microarray analysis, I'm not sure I have grasped the major differences in how they compute differential expression.

Also, how does this manifest in the gene calls? Is one more stringent?

Answer: This is not a comprehensive answer but, in a nutshell, I could say the following about each:

**Limma (as it was/is used for microarrays):**

After background correction, quantile normalisation, and log base 2 transformation, limma fits a linear model to your data and then performs statistical comparisons via linear regression

**DESeq for RNA-seq**

DESeq2 performs a 'geometric normalisation' and then comparisons are generally conducted using the Wald test. DESeq2 also allows you to transform your normalised counts using a regularised log transformation, which deals with transcripts of low counts, or a variance-stabilised transformation, which deals with the high variability of low/high counts that you get from RNA-seq data.
Question: I am trialing discoSNP++ as part of a bacterial GWAS pipeline and am seeking some clarification on the genotypes in the multisample VCFs. Similar to the Pseudomonas example provided in the VCF_creator user guide (pages 4-5), we see a large number heterozygous genotypes despite all the samples being from haploid organisms (e.g. 0/0, 0/1, 1/1). How should we interpret these heterozygous reads from our bacterial sequence data (paired-end data provided in the fof_reads1.txt and
fof_reads2.txt structure as described in Case 4 of the discoSNP user guide). Any guidance appreciated!

Answer: Hi

The genotyping is an option that can be switched off (`-n` option) when working on non diploid species. Computing the genotype on haploid species is meaningless.

However, 0/1 results may warn you (depending on the effective coverage of each allele) as this may reflect the existence of approximate repeats in the genome `-----A------//------T-------` that may be seen as SNP variants while they are not.

1/1 results are expected as, with no reference genome, the "reference" allele is randomly chosen and an homozygous variant may fall in the other allele.

Hope this helps;
Pierre
Question: I would like to know if there are tools or methods for adding a specific sequence into a fastA file at a specific position.

I would like to modify the genome, by adding a specific sequence in a specific location on the genome.

thanks

Assa

Answer: I did [this script][1] for one of my project some years ago. It allows you to look at bases at specific positions but also to modify your genome using another fasta file for the changes.

python lookmod_genome.py -m modify -f file.fasta -c changes_file.fasta

The `changes_file.fasta` should contains all changes you want to make using 4 key words (insertion, deletion, add and remove) and should looks like this :

>insertion:chr1:25:26
GCTAGCTAGC
>deletion:chr4:40:50
>add:your_chromosome_name
GTCGATCGTCATGGTT
>remove:your_chromosome_name

I used it quite a lot but it's self made so check the result with the `look` mode

**Use Python 2 or modify the script to be Python 3 resilient**

[1]: https://github.com/Baboon61/LAM/blob/master/lookmod_genome.py
Question: Hello,

I need to calculate the z-scores from meta-GWAS summary stats. The association file has the following header

rs_id Allele1 Allele2 Freq1 Effect StdErr P.value

The 'effect' is the Beta coefficent. My understanding is

z-score=Effect/StdErr

Is that correct? Also, are there any R packages that can do this?

Many thanks

Answer: You are correct. This can be very easily done with basic AWK or R:

AWK:

awk '{if(NR == 1) print \$1; else if(NR > 1) print $0, \$5/\$6}' infile > outfile

R:

```
Infile<-read.table("infile.txt", head=T)
infile$z-score<-infile$Effect/infile$StdErr
head(Infile, 1)
```

rs_id Allele1 Allele2 Freq1 Effect StdErr P.value z-score
Question: Dear all,

during a step of a CNV analysis on cancer genomes, I would like to be able to COUNT and DISPLAY along the CHROMOSOME AXIS the NUMBER of READS in 10KB REGIONS of GERMLINE and CANCER GENOMES.

I would like to ask you please for your suggestions about :

-- any fast function/algorithm that COUNTS the reads from a BAM file in 10kb windows ?

-- a package/function that DISPLAY the COUNTS along the chromosome axes ?

many thanks,

bogdan

Answer: You should be able to do it with `bedtools`. First generate a BED file of 10kb regions with `bedtools makewindows`. Then calculate coverage for those regions with `bedtools coverage`.

Alternatively, you can try `deepTools bamCoverage` to generate a coverage track (you can specify the size of bins).

Of course, whatever package you are using for CNV calling should have an option to generate some sort of graphical output.
Question: Hi
I have a very simple question.
Is it possible to run a python script in a snakefile, without specifying an output/input?
I have a python script that renames some of my generated output. However, in the script is already specified input files.
Is there a way to just run the script?
My thought we be:

rule pythonscript:
shell:
"script.py"

Is this possible?

Answer: Not knowing the exact input or output at code time is a relatively common thing in bioinformatics, and I think snakemake handles it rather well.

Let's assume `taxonomy.txt` is formatted as below.

[~/Data/scratch/tmp/biostar/checkpoint]$ cat taxonomy.txt
bin	unique	multi	tax
90-20-09-2018.001	25	15	Lactobacillus
90-20-09-2018.003	24	0	Streptococcus
90-20-09-2018.002	15	0	Lactobacillus_2

There are many ways to accomplish it, and below is an example usage of `checkpoint` (https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#data-dependent-conditional-execution ).

[~/Data/scratch/tmp/biostar/checkpoint]$ cat example.py
import os, csv

def aggregate(wildcards):
"""aggregate paths from `read_input_txt` and return as list"""
with open(checkpoints.read_input_txt.get(**wildcards).output[0], 'r') as fin:
return [loc.rstrip() for loc in fin]

rule:
input: aggregate

checkpoint read_input_txt:
"""implement logic to turn `taxonomy.txt` into a list of desired files to create"""
input: 'taxonomy.txt'
output: 'files_to_be_created.txt'
run:
with open(input[0], 'r') as fin, open(output[0], 'w') as fout:
reader = csv.DictReader(fin, delimiter='\t')
writer = csv.writer(fout)
for data in reader:
writer.writerow([os.path.join('binned', data['bin'].split('.')[0], data['tax'])])

rule create_file:
"""implement logic to create the actual file"""
output: touch('{prefix}/{file}')

After `snakemake -s example.py`,

[~/Data/scratch/tmp/biostar/checkpoint]$ tree binned/
binned/
└── 90-20-09-2018
├── Lactobacillus
├── Lactobacillus_2
└── Streptococcus

1 directory, 3 files

Hope you'll find this helpful and able to expand the example into a working solution.
Question: Dear Biostars, Hi

I have a **volcano** plot (obtained from edgeR). I also have **some selected annotated genes** that I like to highlight them by **showing** only their name on that plot.

I have used the [valuable script/code][1] from Biostars (thank you @WouterDeCoster and @venu and others).

As most of the lines of the first column in my *counts.matrix* is empty (I have only about 15 names), I received some errors. For example, in the "Head" example below, I have only **PAC1** and **Gh** as Gene name and the other lines have "**space**" instead of any Gene name!

would you please check my script and my counts.matrix example and help me in this regard? Thanks

NOTE:

1- Head of my counts.matrix file:

Gene	external_gene_name	logFC	logCPM	P.Value	FDR
PAC1	TRINITY_c0_g2_i4	-13.8732718193368	4.61855606743036	3.03811154383156e-30	1.87897780840196e-24
TRINITY_c5_g3_i4	13.2697070220115	4.01618042373378	1.21936893870205e-27	3.7707094407506e-22
GH	TRINITY_c9_g2_i2	11.2269211011449	3.82178000481344	3.20153046019011e-26	6.60015780727773e-21
TRINITY_c0_g1_i1	12.855465908739	3.60335598332695	2.98018910927419e-24	4.60788644555924e-19
TRINITY_c0_g1_i1	-11.8727923903381	2.62735890209779	1.81499535997098e-20	2.24503673057178e-15
TRINITY_c6_g1_i6	10.3243429332541	2.89490936855587	2.16106621247948e-19	2.22758743227662e-14
TRINITY_g3_i2	-12.0236467730875	2.77706530918009	2.63166646685828e-19	2.32514875441625e-14
TRINITY_c3_g2_i6	-11.2400737828536	2.00142608642608	1.40238473598696e-18	1.0841643566014e-13
TRINITY_c1_g1_i1	-9.45052158373113	2.158090406858	1.66695714263313e-18	1.14551257449686e-13

2- my code:

# Load packages
library(dplyr)
library(ggplot2)
library(ggrepel)

# Read data
data <- read.table("trans.counts.matrix.J_vs_M.edgeR.DE_results", header=TRUE)

data$threshold = as.factor(data$FDR < 0.001)

data = mutate(data, sig=ifelse(data$FDR<0.001, "FDR<0.001", "Not Sig"))

p = ggplot(data, aes(logFC, -log10(P.Value))) +
geom_point(aes(col=sig)) +
scale_color_manual(values=c("red", "black"))

p
p+geom_text_repel(data=filter(data, FDR<0.001), aes(label=Gene))

[1]: https://www.biostars.org/p/203876/

Answer: Change this line of code

data <- read.table("trans.counts.matrix.J_vs_M.edgeR.DE_results", header=TRUE)

to

data <- read.delim("trans.counts.matrix.J_vs_M.edgeR.DE_results", header=TRUE, na.strings=c("","NA"))

This adds `` to empty cells which will be ignored in the plot and only gene names will be added.
Question: How to get the gene length from a list of Ensembl IDs using biomaRt for instance (or any R-based method without having to download a separate annotation file first)?

Since `gene_length` attribute does not exist in biomaRt, is there a better alternative than using `start_position` and `end_position` attributes, then substracting the 2 values like as follows:

library(biomaRt)
ensembl_list <- c("ENSG00000000003","ENSG00000000419","ENSG00000000457","ENSG00000000460")
human <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
start_pos = getLDS(attributes = "ensembl_gene_id", filters = "ensembl_gene_idl", values = ensembl_list , mart = human, attributesL = "start_position", martL = human, uniqueRows=T)
end_pos = getLDS(attributes = "ensembl_gene_id", filters = "ensembl_gene_idl", values = ensembl_list , mart = human, attributesL = "end_position", martL = human, uniqueRows=T)
gene_L <- merge(start_pos, end_pos, by.x="Gene.stable.ID", by.y="Gene.stable.ID")
gene_L$Length <- gene_L$Gene.end..bp. - gene_L$Gene.start..bp.

Answer: library (EDASeq)
ensembl_list <- c("ENSG00000000003","ENSG00000000419","ENSG00000000457","ENSG00000000460")
getGeneLengthAndGCContent(ensembl_list, "hsa")

Connecting to BioMart ...
Downloading sequences ...
length gc
ENSG00000000003 4535 0.3995590
ENSG00000000419 1207 0.4159072
ENSG00000000457 6883 0.4117391
ENSG00000000460 5967 0.4298643

@OP: I think code is unnecessarily complex (IMO) for the logic in OP and it can be shortened further by:

> library(biomaRt)
> ensembl_list <- c("ENSG00000000003","ENSG00000000419","ENSG00000000457","ENSG00000000460")
> human <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
> gene_coords=getBM(attributes=c("hgnc_symbol","ensembl_gene_id", "start_position","end_position"), filters="ensembl_gene_id", values=ensembl_list, mart=human)
> gene_coords$size=gene_coords$end_position - gene_coords$start_position
> gene_coords

hgnc_symbol ensembl_gene_id start_position end_position size
1 TSPAN6 ENSG00000000003 100627109 100639991 12882
2 DPM1 ENSG00000000419 50934867 50958555 23688
3 SCYL3 ENSG00000000457 169849631 169894267 44636
4 C1orf112 ENSG00000000460 169662007 169854080 192073

However these number are not matching with those from EDASeq. I checked the coordinates of gene TSPAN6 ( ENSG00000000003) on NCBI. In NCBI gene (7105) start and stop respectively are: 100627108 and 100637104. Start coordinate on NCBI is off by 1 base from biomart results and stop coordinate is farther away (on chromosome x) compared to end position on NCBI. When viewed (the same gene) in IGV (hg38), gene start and stop positions are: chrX:100,625,108-100,639,104. In general, for a given gene, IGV display coordinates are 2kb up and down stream of gene start and end. Thus hg38 coordinates (from IGV) and GRCh 38 match (quick way). EDAseq probably is doing sum of exon lengths and calculating GC content based on that.
Question: Hi,

I have a data.frame of 17.7M rows, which I would like to plot as a heatmap (heatmap.2).

I would first like to create a table to see how many elements I have between (0,1], (1,2], (2,3]... ,(9,10]. I know that I have ~3500 elements without 'NA' and ~ 700 below 1. How can I do a complete overview of the data? Can I do it with `table`?

I would than like to bin the same data into specific (factorial?) values, so that I can use my own colors legend for each range.

This is my data structure:

```
> head(filtered1000)

start start.1 Value
1 876001 1403001 9.910803
2 3079001 3081001 9.891197
3 834001 836001 9.813543
4 1239001 1241001 9.794319
5 2777001 2780001 9.775171
6 1712001 1714001 9.727626
```

I know it is probably possible to manually adjust the `cut` command into the right breaks, but is there a more "subtle", automated way of calculating the number of each elements in each category and apply them to specific bins so that I will be able to set a certain color to a specific bin?

```
> breaks <- c(0,1,2,3,4,5,6,7,8,9,10)
> bins<-cut( filtered1000[,3], breaks )
summary(bins)

(0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] NA's
2798 366 172 26 45 5 11 18 25 31 17762728
```

How can I than incorporate these breaks into my data.frame, so that I can plot (ggplot2) it with specific colors?

Thanks

Assa

Answer: It would be simpler to use `findInterval()`.
Question: I just got my RNA seq data.

CZ1_S6_L001_R1_001.fastq

CZ1_S6_L002_R1_001.fastq

CZ1_S6_L003_R1_001.fastq

CZ1_S6_L004_R1_001.fastq

I only know there are from 4 different lanes during sequencing.

What are they and do I need to combine these fastq files into one? And how?

Thanks~~~

Answer: You first might want to QC your data lane by lane to compare lane effects.

Based on the filenames I get the idea your data is single ended? In that case you can just concatenate the files for downstream analysis.

e.g.

cat CZ1_S6_L001_R1_001.fastq CZ1_S6_L002_R1_001.fastq CZ1_S6_L003_R1_001.fastq CZ1_S6_L004_R1_001.fastq > CZ1_S6_merged_R1_001.fastq
Question: DESeq2 vignette states

> The point of these two transformations, the VST and the rlog, is to
> remove the dependence of the variance on the mean, particularly the
> high variance of the logarithm of count data when the mean is low.

and the documentation of rlog explains

> The transformation is useful when checking for outliers or as input for machine learning techniques such as clustering or linear discriminant analysis

I understand that "checking for outliers" means checking for outliers via a PCA plot (or something similar).

Why is minimizing differences (between samples) for rows with low counts important for the PCA plot?

Why does the variance have to be independent of the mean (homoscedasticity) for that?

Answer: In variance based analyses, like PCA, clustering and LDA, the results are driven by those features with the highest variance. In count-based data, like RNA-seq, there is a relationship between the mean and the variance - higher mean = higher variance when the data is on the linear scale.If you were to run something like PCA on the linear scale, you would simply find that the result was dominated by the random noise in the high mean features.

However, because of the discrete nature of the data, on the log scale, the variance is higher in the low count features: 2 reads is twice as many as 1, and 1 read is infinitely more reads than 0, but 101 reads is only 1% more than 100. Thus, without some kind regularization, your PCA will be dominated by very small changes in low count features.
Question: I have some older VCF files that don't have the contig length set in the VCF header. This means that Picard and some other tools that are very strict with the VCF spec won't accept them.

The contig entries in the header should be

##contig=
##contig=

but they are

##contig=
##contig=

I know that I can manually fix this by doing the following steps.

- unzipping the file
- extracting the header
- lookup the contig lenghts in a fasta.fa.fai file
- adding the lenght to the contigs records in the header with vim
- re-header with bcftools
- bgzip and tabix the the re-headered vcf file

In my hands this works but if you make a slight VIM copy paste error you will have spend a lot of time reheadering and bgzipping a large VCF file for nothing.

Therefore I would like to have a more robust automatic solution where I just give the VCF file and the reference genome file and the header is automatically fixed and a new bgzipped VCF file is written out.

Is there a tool that can add the contig lenghts to the VCF header and write out a new bgzipped VCF file?

Answer: A select all variants with the gatk framework also adds the contig lenght to the contigs in the header. Upside is that it's a one liner with a off the shelf tool:

`gatk-framework -T SelectVariants -V input.vcf.gz -R /ref.fa -o input_header_fixed.vcf.gz`

Downside side is that the full VCF is inflated to the GATK Variant/GenotypeContext objects and rewritten to a new VCF. So it's slower then AWK and or a c/c++ tool solution, and the new VCF is slightly different (different order or precision of Variant/Genotype attributes).
Question: Hi!

I was wondering whether it is feasible to work with multiple types of data from same patient in TCGA to reveal mechanism of a cancer in that particular patient using tumor and adjacent healthy tissue controls. I am planning to analyse each data type individually and merge findings in a higher-level (e.g. pathway most probably).

I have been reading on it a lot but could be able to find satisfactory answers to following:

1. On TCGA data quality, what do you think on the reliability of data for individual analysis? I have received different opinions on this matter, so I'm a bit confused right now.
2. For all types of data, since tumor and adjacent tissue controls are in the same batch in TCGA, how should I check for other possible biases? I guess batch effect correction becomes irrelevant in this case. This question particularly becomes relevant in array based methods, such as HumanMethylation450k. If tumor and control scanned in different arrays, how can I be sure that there is no additional bias coming from array scanning (except probably expecting low DMR count)?

I really appreciate your valuable input in these matters.

Thanks a lot!

Answer: On point 1, there is not really a general answer to your question. If your analysis demands high-quality data (and not all do), then doing some QC yourself is definitely necessary. On point 2, not all data types for all tumor types include "normal" controls, so the points you bring up may not be relevant. As for technical variation, that is something largely out of your control once the data are collected. Using formal hypothesis testing with multiple samples is really one of the best ways to approach datasets that can deal with variation in the data in an efficient manner (statistically speaking).