cannin/biostars_qa huggingface dataset:
Question: I was searching for a function in R to perform a chi-squared test for hardy-weinberg equilibrium.

Strait-forward would be to compute a chi-squared-value by my own and then take a look at a p-value with a given degree of freedom:

```
nAA=18
nAa=104
naa=109
n=nAA+nAa+naa
pA=(2*nAA+nAa)/(2*n)
pa=1-pA
obs= c(nAA,nAa,naa)
exp=c(n*pA*pA,2*n*pA*pa,n*pa^2)
sum(((obs-exp)^2)/exp)
[1] 1.001318
1-pchisq(1.001318,df=1)
[1] 0.3169918
```

My degrees of freedom are 1, as number of genotypes = 3 and of alleles = 2 (df=3-2=1). So, p-value is 0.31699, I am happy, everything is in HWE. Then I decided to have a look if there is a function in R for that within standard packages and have found the following:

```
chisq.test(obs, p = exp/sum(exp))
Chi-squared test for given probabilities data: obs X-squared = 1.0013, df = 2, p-value = 0.6061
```

It seems that `chisq.test` is not suitable for HWE as it computes the df as 2. Are there any functions that can compute the `chisq.test` for HWE problems?

Thank you.

Answer: A quick search turns up the HardyWeinberg R package. Specifically, HWChisq() performs the test you want.
Question: Hi

I have a compiled segmentation file (nocnv) for almost 1000 samples. I want to calculate the copy number aberrations for all the samples.

I tried using GISTIC, but I believe it provides only the consensus copy number variations. Is there a better way to calculate copy number aberration from segmentation file.

Here are few lines from the segmentation file I have:

```
Sample Chromosome Start End Num_Probes Segment_Mean
TCGA-E2-A1IH-01A-11D-A13J-01 1 3208470 42629316 20791 -0.2864
TCGA-E2-A1IH-01A-11D-A13J-01 1 42631290 42692913 60 0.0183
TCGA-E2-A1IH-01A-11D-A13J-01 1 42701180 42918495 113 0.3411
TCGA-E2-A1IH-01A-11D-A13J-01 1 42925163 43153569 131 -0.0548
TCGA-E2-A1IH-01A-11D-A13J-01 1 43153652 43599777 288 -0.3653
TCGA-E2-A1IH-01A-11D-A13J-01 1 43600151 44642193 575 -0.0401
```

Thanks

Meenakshi

Answer: You'll want to define some threshold for `|Segment_Mean| > X` because I would probably not consider all segments to be true copy number alterations. Filtering based upon length and number of probes is probably also a good idea.

I don't have a precise recommendation of what threshold to use, but I would try plotting density / histogram values for the Segment_Mean to guide your decisions.
Question: Hi all,

I have a list of more than 100 SNPs (rsXXXXXX) and I would like to obtain the different allele frequencies that each of them shows in each of the 1000Genomes populations (if possible not manually...). Is there any tool, R package etc, that allows to download them all at once?
I've thought that I could maybe obtain those frequencies from the UCSC database or directly from the 1000Genomes database, but I'm open to any suggestions.
Thank you very much in advanced!

Answer: Solution 1: The raw variant call data can be downloaded from http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ . Once you have those files, the INFO column in each VCF file contains superpopulation allele frequencies, and there are a bunch of tools which can look up the INFO column entry for a particular rsID. The one complication is that, since there's a separate VCF file per chromosome, you may first need to figure out which rsIDs are on which chromosomes.

Solution 2: With plink 2.0, http://www.cog-genomics.org/plink/2.0/resources#1kg_phase3 provides a single merged dataset containing all chromosomes (download the boldfaced links, then rename phase3_corrected.psam to all_phase3.psam). Then,

plink2 --pfile all_phase3 vzs --extract [your list of rsIDs] --export vcf

can then be used to export a VCF with only the rsIDs you care about; the precomputed superpopulation allele frequencies will be in the INFO column of this freshly generated VCF. You can also define your own populations with --keep and compute allele frequencies on the fly with --freq.
Question: Hi all,

I have a FASTA file which contains protein sequences of a load of genes from D. melanogaster, and need to split the file into multiple FASTAs, one gene per file. What's the best way to go about this? Ideally each file will be named with the name of the gene (UniProt ID).

Example, the sequence below wants to be split into 2 files, the first called O46197.fasta, the second called Q9VUQ5.fasta

Many thanks!

>sp|O46197|A29AB_DROME Accessory gland protein Acp29AB OS=Drosophila melanogaster GN=Acp29AB PE=2 SV=2
MYASNLLYLLALWNLWDLSGGQQDIPNGKATLPSPQTPQNTIDQIGINQNYWFTYNALKQ
NETLAIIDTMEMRIASSLLEFKAQMEIQLQPLKIIMRHHASNIKASNNIKMRRFEKVGSR
HFHIEKNLMQTWFEAYVTCRKMNGHLANIQDEMELDGILALAPNNSYWIDISKLVENGGT
FVSTLTGREPFFVKWKSNQDTKKKNQCVYIYAKEMSYDECFEKKSFVCQADQWA

.
>sp|Q9VUQ5|AGO2_DROME Protein argonaute-2 OS=Drosophila melanogaster GN=AGO2 PE=1 SV=3
MGKKDKNKKGGQDSAAAPQPQQQQKQQQQRQQQPQQLQQPQQLQQPQQLQQPQQQQQQQP
HQQQQQSSRQQPSTSSGGSRASGFQQGGQQQKSQDAEGWTAQKKQGKQQVQGWTKQGQQG
GHQQGRQGQDGGYQQRPPGQQQGGHQQGRQGQEGGYQQRPPGQQQGGHQQGRQGQEGGYQ
QRPSGQQQGGHQQGRQGQEGGYQQRPPGQQQGGHQQGRQGQEGGYQQRPSGQQQGGHQQG
RQGQEGGYQQRPPGQQQGGHQQGRQGQEGGYQQRPPGQQQGGHEQGRQGQEGGYQQRPSG
QQQGGHQQGRQGQEGGYQQRPSGQQQGGHQQGRQGQEGGYQQRPSGQQQGGHQQGRQGQE
GGYQQRPPGQQPNQTQSQGQYQSRGPPQQQQAAPLPLPPQPAGSIKRGTIGKPGQVGINY
LDLDLSKMPSVAYHYDVKIMPERPKKFYRQAFEQFRVDQLGGAVLAYDGKASCYSVDKLP
LNSQNPEVTVTDRNGRTLRYTIEIKETGDSTIDLKSLTTYMNDRIFDKPMRAMQCVEVVL
ASPCHNKAIRVGRSFFKMSDPNNRHELDDGYEALVGLYQAFMLGDRPFLNVDISHKSFPI
SMPMIEYLERFSLKAKINNTTNLDYSRRFLEPFLRGINVVYTPPQSFQSAPRVYRVNGLS
RAPASSETFEHDGKKVTIASYFHSRNYPLKFPQLHCLNVGSSIKSILLPIELCSIEEGQA
LNRKDGATQVANMIKYAATSTNVRKRKIMNLLQYFQHNLDPTISRFGIRIANDFIVVSTR
VLSPPQVEYHSKRFTMVKNGSWRMDGMKFLEPKPKAHKCAVLYCDPRSGRKMNYTQLNDF
GNLIISQGKAVNISLDSDVTYRPFTDDERSLDTIFADLKRSQHDLAIVIIPQFRISYDTI
KQKAELQHGILTQCIKQFTVERKCNNQTIGNILLKINSKLNGINHKIKDDPRLPMMKNTM
YIGADVTHPSPDQREIPSVVGVAASHDPYGASYNMQYRLQRGALEEIEDMFSITLEHLRV
YKEYRNAYPDHIIYYRDGVSDGQFPKIKNEELRCIKQACDKVGCKPKICCVIVVKRHHTR
FFPSGDVTTSNKFNNVDPGTVVDRTIVHPNEMQFFMVSHQAIQGTAKPTRYNVIENTGNL
DIDLLQQLTYNLCHMFPRCNRSVSYPAPAYLAHLVAARGRVYLTGTNRFLDLKKEYAKRT
IVPEFMKKNPMYFV

Answer: try the following awk script:

awk -F '|' '/^>/ {F=sprintf("%s.fasta",$2); print > F;next;} {print >> F;}' < yourfile.fa
Question: Hello!

I need to make contigs from the Ion Torrent data. This is a microbial genome fragment library. It supposed to be a set of 400 bp read. Due to quality reasons (I used FastQC to check) the reads are kind of good for the first 170-200bp. Our core sequencing facility releases the data as *.FASTQ files and the adaptors are supposed to be trimmed.
I am not very experienced with NGS data analysis. The only program I used to clean my date was Trimmomatic but at that time I worked with Illumina reads.
Could you suggest me something like Trimmomatic for my Ion-torrent data, please?

Thanks for your tips.

Answer: Optimal would be to use the ubam file straight from the sequencing machine with Spades assembler. You lose valuable information (that Spades utilizes) when you convert the ubam files to fastq.
Question: Hi all, was wondering what does it mean by secondary alignment in BAM flag?

What is been defined as unmapped reads?

**Option used**:

samtools view -F 260 bam_file

Thanks!

Answer: Unmapped reads are in the BAM file but have no valid assigned position (N.B., they may have an assigned position, but it should be ignored). It's typically the case that a number of reads can't be aligned, due to things like sequencing errors, imperfect matches between the DNA sequenced and the reference, random e. coli or other contamination, etc..

A secondary alignment occurs when a given read could align reasonably well to more than one place. One of the possible reported alignments is termed "primary" and the others will be marked as "secondary".
Question: This is my bed file for all the exon coordinates ..small subset so i want to take out all the exons that of a given gene let say i have gene in chr 1 which starts from `chr1	11868	12227` so i want to parse out all the exons that comes in between `11868	12227`

this is my small subset

cat exon.bed | head -10
chr1	11868	12227	+	exon
chr1	11871	12227	+	exon
chr1	11873	12227	+	exon
chr1	12009	12057	+	exon
chr1	12178	12227	+	exon
chr1	12594	12721	+	exon
chr1	12612	12697	+	exon
chr1	12612	12721	+	exon
chr1	12612	12721	+	exon
chr1	12974	13052	+	exon

How do i parse out , i use mostly R and bit of shell script but I m not sure if i can use R , may a few lines of perl or shell script can help me solve my problem.

Any help or suggestion would be highly appreciated

Answer: Via BEDOPS `bedops -n` and Unix I/O streams:

$ echo -e "chr1\t11868\t12227" | bedops -n 1 exon.bed - > answer.bed

Or, if you have your genes in a BED file called `genes.bed`:

$ bedops -n 1 exon.bed genes.bed > answer.bed

If you have your genes in some other format, like GFF or GTF, you can use `gff2bed` or `gtf2bed`, e.g.:

$ bedops -n 1 exon.bed <(gff2bed < genes.gff) > answer.bed

Or:

$ bedops -n 1 exon.bed <(gtf2bed < genes.gtf) > answer.bed

The file `answer.bed` will contain exons that do not overlap a gene annotation.
Question: First, I don't understand what it means by partial correlation matrix. I know correlation matrix, though.

Second, assuming partial correlation matrix is similar to correlation matrix, then how can you refer the connections (edges) between the dots (genes) and how to resolve the directionality?

Third, I am using the R package "parcor". I have read the R manual, some reference papers, but still don't know how to infer a network from the package. Can somebody give a simple tutorial?

Thanks much in advance.

Answer: The partial correlation between X and Y is the correlation between the residuals of X and Y after linear regression to eliminate the effect of confounding factors. For example, assume you have two independent variables W and X and a dependent variable Y. To compute the partial correlation between X and Y, you regress X on W to get the residual rX and regress Y on W to get residual rY. The partial correlation is then the correlation between rX and rY. It allows to answer questions like how much X contributes to Y when the effect of W is removed.
As the partial correlation can be viewed as a similarity measure, a partial correlation matrix can be seen as the adjacency matrix of a graph where the partial correlations represent edge weights. Since it's a correlation matrix, it is symmetric. A symmetric adjacency matrix implies an undirected graph so you need additional information to build a directed graph out of it.
As for the practical aspect, once you have your matrix, you can use the [graph_from_adjacency_matrix()][1] function from the [igraph package][2] to get the coresponding graph structure.

[1]: http://igraph.org/r/doc/graph_from_adjacency_matrix.html
[2]: http://igraph.org/r/
Question: In biopython `feature.location.end` can be equal to `next_feature.location.start`. For example:

```
type: TRANSMEM
location: [187:208]
qualifiers:
Key: description, Value: Helical. {ECO:0000255}.
type: TOPO_DOM
location: [208:411]
qualifiers:
Key: description, Value: Extracellular. {ECO:0000255}.
```

Although there is some biological ambiguity over this example (residue 208), in others there is not. Hence I ask which domain do the residues that are overlapped truly belong to?

Answer: Devon is right that a protein might well be annotated with overlapping domains, however in this case the domains in your example do NOT overlap. Biopython uses Python style slicing notation, so [187:208] and [208:411] do NOT overlap. e.g.

```
>>> example = "0123456789"
>>> example[3:6]
'345'
>>> example[6:9]
'678'
```

Also beware that Biopython and Python use zero-based counting, rather than the one-based counting you may be more used to. Note SwissProt/UniProt annotation files use one-based counting in their plain text and XML file formats.
Question: Good afternoon :)

This is more of a theory question. My general pipeline is BBDuk.sh -> STAR align -> salmon quant. I am using 3' RNA-seq (Lexogen QuantSeq FWD prep).

During BBDuk.sh, this does a poly a tail trim, as well as a adapter trim (all happy here).

I then align the trimmed reads to the reference genome using STAR (again, all happy here).

This is where some differences occur with Salmon versions.

Using salmon v0.10.0, i get expected outputs and quantified files.

/nethome/bdy8/programs/salmon-0.10.0_linux_x86_64/bin/salmon \
quant \
-t /nethome/bdy8/apal_genome/version3.0/Apalm_gffread_for_salmon.fasta \
-l SF \
-a /scratch/projects/transcriptomics/ben_young/DHE/tagseq/host/aligned/'"$PALPAL"'/'"$PALPAL"'_Aligned.toTranscriptome.out.bam \
-o /scratch/projects/transcriptomics/ben_young/DHE/tagseq/host/salmon/'"${PALPAL}"'_salmon

However, using salmon 1.40 I experience thew following error message

/nethome/bdy8/programs/salmon-latest_linux_x86_64/bin/salmon quant \
-t /nethome/bdy8/apal_genome/Apalm_gffread_for_salmon.fasta \
-l SF \
--noLengthCorrection \
-a /scratch/projects/transcriptomics/ben_young/POR/tagseq/host/aligned_reads/bagnumber-apal-1009/bagnumber-apal-1009_Aligned.toTranscriptome.out.bam \
-o /scratch/projects/transcriptomics/ben_young/POR/tagseq/host/test

> SAM file says target evm.model.Sc0a5M3_402_HRSCAF_756.4.1.5f5b2bc4 has
> length 536, but the FASTA file contains a sequence of length [538 or
> 537]

What I am thinking is happening is that this is resulting due to the poly a tail trimming in the BBDuk.sh stage. Apart from this I do not know what is happening and was wondering if anyone else has any ideas on this or a fix to ignore this mismatch of 1.

If any more information is needed please let me know and I will supply it.

Answer: To provide an answer to this query, it seems STAR is doing something during the indexing step which is causing a slight mismatch for 23 of the transcripts. If you read above replies to the original post you will see the workflow that was done to identify the problem transcripts and the various lengths of them in different programs :).

I will be opening a question with Alex (STAR author) about identifying why this may be happening.
Question: I am using SNPRelate for PCA analysis. Its using default color for different populations but I want to color them according to me. Plotting commands are like this:

```
> plot(tab$EV2, tab$EV1, col=as.integer(tab$pop),cex=1.2,pch=20,
+ xlab="eigenvector 2", ylab="eigenvector 1")
> legend("topleft", legend=levels(tab$pop), cex=1,pch=20, col=1:nlevels(tab$pop))
```

Head of the input file is like this:

```
sample.id pop EV1 EV2
1 A1 POP_I -0.10172849 0.03619405
2 A2 POP_I -0.15951814 0.08234857
3 A3 POP_I -0.15632495 0.08180843
4 A4 POP_I -0.09679447 0.07981108
5 A5 POP_I 0.11362360 -0.03186038
6 A6 POP_I 0.05594095 -0.05498351
```

Any idea about it?

Answer: I think you want to use R's [palette][1] command.

[1]: https://stat.ethz.ch/R-manual/R-devel/library/grDevices/html/palette.html
Question: I have used this piece of script to draw the PCA based on eigenvalue and eigenvector

percentage <- round((eigenval/(sum(eigenval))*100), 2)
percentage <- as.matrix(percentage)
percentage <- paste0(names(percentage), " (", percentage, "%)")
Names <- c ("mn27hd", "mdkk987", "mnsdnu83", "sjednu83", "bjeo972s")
pop.colour <- c("blue", "red", "green", "orange", "brown")
ggplot(eigenvec, aes(x=PC1, y=PC2, colour=pop.colour, label=Names)) +
geom_point(size=3) + geom_text(aes(label=Names), hjust=0, vjust=-0.5) + xlab(percentage[1]) + ylab(percentage[2])

![enter image description here][1]

[1]: /media/images/c26677bb-7e25-433a-a8d7-eb7b43ad

how can I:

1. change the xlab and ylab names to show as: PC1 (36.61%), PC2 (33.7%)?

2. move only the pink label slightly left-side cause when I do it generally then the two left labels go outside the image or too close to the edge?

3. change the legend name and categories from the colour to the population name?

Answer: I like to use `labs` instead of `xlab` and `ylab` and add the information to the data.frame directly. Do you only have 4 data point and you are certain that the name are corresponding to the correct PC value? If so do
```
percentage <- round((eigenval/(sum(eigenval))*100), 2)
percentage <- as.matrix(percentage)
percentage <- paste0(names(percentage), " (", percentage, "%)")
eigenvec$Pop <- c ("mn27hd", "mdkk987", "mnsdnu83", "sjednu83", "bjeo972s")

pop.colour <- c("blue", "red", "green", "orange", "brown")
ggplot(eigenvec, aes(x=PC1, y=PC2, colour=Pop, label=Pop)) +
geom_point(size=3) + geom_text(hjust=-0.1, vjust=-0.5) +
scale_color_manual(values = pop.colour) +
labs(x=paste0("PC1 (",percentage[1],")"), y= paste0("PC2 (",percentage[2], ")"), color = "Population")

```

I don't know have mock data to play with so there might be bug to the code. You can modify `hjust` in `geom_text` to adjust the horizontal shift
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: command line:

```
egrep "^[0-9XY]" file
```
EDIT: originally, this used the same regex as for the R-based example (`^[0-9XY]$`). This won't work because the full line of the text file contains more characters (such as the coordinates...). Thanks to Alex for pointing this out.

R (because this post is tagged with it):

```
# assuming your data frame with the coordinates looks like this
df <- data.frame(chr = c("1","2","X","Y", "GL000220.1"),
start = c(1,20,30,40,50),
end = c(11, 21, 31, 41, 51)
)

subset(df, grepl("^[0-9XY]$", chr))
```
Question: Hi Everyone

I converted the gtf file I got from stringtie to bed.
The bed still contains the stringtie gene and transcripts code (same thing if I used cufflinks instead) , is there a way to post-process this bed file to have real gene names and transcripts not stringtie nor cufflinks ids:

chrY	23451597	23451912	STRG.80985.1	0	.	23451597	23451912	0	1	315,	0,
chrY	23460655	23460996	STRG.80986.1	0	.	23460655	23460996	0	1	341,	0,
chrY	23464508	23464868	STRG.80987.1	0	.	23464508	23464868	0	1	360,	0,

Thanks

Answer: You could use BEDOPS `bedmap` to annotate regions with genes (or their names):

$ bedmap --echo --echo-map-id regions.bed genes.bed > answer.bed

If you edit your question to note the genome you're working with, I can suggest where to get annotations to build the `genes.bed` file.
Question: Hi I have a question, I have Drosophila genome spike-in in a chip-seq experiment. Does anybody have experience on the best way to normalised and process those samples to generate bw and count matrix? ( I am using deeptools)

I am following the instruction of the vendor:

map reads to the sample genome (human, mouse, etc.)
map reads to the DM genome
count unique DM sequence tags. Identify sample with least amount of tags
Count DM tags from other samples. Create a normalization factor = lowest_sample/sample_of_interest

then I use the normalisation factor ${SPIKE} in deeptools as follow:

bamCoverage -b ${SAMPLE}/${SAMPLE}.nodup.norm.bam -o /scratch/${PBS_JOBID}/${SAMPLE}/${SAMPLE}.bw \
-p 1 --normalizeTo1x ${genomesize} \
--scaleFactor ${SPIKE} \
--ignoreForNormalization chrX chrY \
-bl ${ENCODEBED} \
--smoothLength 40 \
--binSize 10 \
--centerReads --extendReads 150

Since I am getting a modest increase, and the vendor results are a modest decrease of signal I am wondering if I am doing something wrong with bamCoverage --scaleFactor argument or if the discrepancy is due to the different processing of the bw between my pipeline (deeptools based) and the vendor one (custom made).

Thanks a lot for your support and advice

Answer: Did the vendor instructions mention that they also normalized to 1x? That could be the problem...

It's possible that you should only be applying the scaleFactor alone and not also normalizing to 1x. If you apply only the scaleFactor (i.e. run that command without the --normalizeto1x argument), would you see a modest decrease of signal as the vendor does?

The [deepTools google group][1] is another great resource for questions like this one, and they may be able to help more.

[1]: https://groups.google.com/forum/#!forum/deeptools
Question: Hello,

If I have TPM values for gene A from the UCSC database which contains multiple isoforms, and lets say I want just 1 output for gene A, can I average or add the TPM values?

If I can't do that is it possible to go back to the raw reads and add them at the raw read level?

I just want to condense genes for my dataset.

Answer: If they are calculated correctly (e.g. with RSEM or similar), the quantification should be divided among the isoforms in a mutually exclusive way.

So if you have

Isoform1 3 TPM - three transcripts
Isoform2 1 TPM - one transcript

The total gene should have 4 transcripts per million (i.e. 4 independent transcripts, 3 from Isoform1 and one from Isoform2.)
Question: Hi all,

I have a bacterial genome that I know has an approximate size of 7.3m. I have both short-reads (Illumina) and long-reads (Pacbio). Sequencing showed good coverage.
I tried to assemble these data using SPAdes, Flye or Canu but I can not obtain a circular chromosome contig. I know there is plasmid contigs in this genome but the best assembly I got was with Flye (and rounds of Pilon) with a total of 8 contigs.

Do you know another assembler tool that I can use? Do you have any advices on how should I perform this assembly?

Thank you very much for your help,

Answer: I'd suggest Unicycler. From my experience it sometimes helps to assemble the long reads first using Flye and feed that as a pre-assembly back into Unicycler using both short and long reads.
Tricycler would also be an interesting option with various rounds of short read polishing afterwards - but no experience so far with that.
Question: Hello.

I was working on a pipeline for genome assembly and I have used the human paired-end NGS data.

I want to compare my denovo assembly to the available human genome which is based on Sanger sequencing mainly.
How can I compare my denovo assembly to the sanger sequencing inorder to check -

How much of the genome was covered in the denovo assembly?

How many regions of the genome have I missed in the denovo contigs?

Do I compare the contigs or scaffolds to the reference which is ordered as chromosomes?

I know Quast is used to test different genome assembly qualities, but what about comparing a denovo assembly to a reference genome???

Answer: Quast also works with a reference genome, just set the reference with -R. According to their manual, Quast has tests that will give you answers to points 1 & 2. As for point 3, Quast has a --scaffolds option which, if given scaffolds will break them back into contigs. You could run it with and without the --scaffolds option to see the difference in the metrics. There is also ALE for looking at assembly metrics although I'm not very familiar with it