cannin/biostars_qa huggingface dataset:
Question: After I got the filtered vcf file of snps with gatk pipeline, I tried to run PCA with gcta.
afterwords, I tried to find the explained variation percentage for each PCs(principal components).
what I'm confused is that the eigenvalue results I get are always the same number as the sample size.
I thought that the amount of possible PCs were the same as the variable dimension I have (which is the number of snps called, about 17000000 in this case) and there should be the same amount of eigenvalues related to it. Of course there could be same eigenvalues and PCs that explain small variation are useless, but isn't it possible that there could be 17000000 PCs and eigenvalues? So I thought when I want to get the explained variation percentage of PC1, I had to divide it by the some of 17000000 eigenvalues. could someone explain why this is wrong and why I always get the same number of eigenvalues as the sample size?

Answer: The principal components are the eigenvectors of the covariance matrix (or correlation, if you scale) of your data. Now, you have two choices: do you want to focus on the covariance between your N samples, or the covariance between your M SNPs? GCTA focuses on the former, whereas you thought it was the latter. GCTA's covariance matrix (aka GRM) is NxN, hence you will get N eigenvectors. You are not reducing the "SNP dimensions", you are reducing the "sample dimensions".
Question: I'm trying the tuxedo protocol as well as analysis using DESeq2, I wanted to see if I can see in the difference of DE genes. DESeq2 uses something called count matrix ,which I generated using both ht-seq and featurecounts .Now I want to know how do I transform those counts to make it usable for DESeq2 analysis

I see the count table from ht-seq and featurecount whats the difference between both the output.? for ht-seq I see a long list of ensemble ID with the counts where as for featurecount there is lot more information .

I have two 4 samples Control vs Test with replicate of each of them ,so I have the all together 4 count table.

So how can i use those counts into my DESeq2 ?

Any help and suggestion would be highly appreciated.

Answer: The output of both HT-seq and featurecounts can *almost* be fed directly in DESeq2.

- **for FeatureCounts :**

Output looks like :

Geneid	Chr	Start	End	Strand	Length	cond1	cond2
SPBC460.05	I	16470	18062	+	1593	1	12	24
SPBC460.02c	II;II	8856;9651	9365;9803	-;-	663	329
SPAC212.11	I	1	5662	-	5662	0	0	0

First you need to remove 5 firsts columns and name the rows as gene names

cond1	cond2
SPBC460.05	12	24
SPBC460.02c	663	329
SPAC212.11	0	0

then just use the `DESeqDataSetFromMatrix` function of DESeq2 to convert it.

- **for HTseqCounts :** see comments below
Question: Hi,
I would like to calculate the percent-identity from a CIGAR string from a BAM/SAM file containing alignments. I want to calculate the PID only for the aligned region, ignoring clipped ends ("H","S"). I can parse the CIGAR in R and get the sums for each letter in the CIGAR, so it's just a conceptual question about the definition and whether or not it is correct:
so given the CIGAR contains characters: "M" (match),"N" (skip),"D" (Deletion),"I"(Insertion), "S" (soft-clip),"H" (hard-clip), of which I ignore "S", "H", and "M" is the total sum of M, N the total of N, etc.:

e.g.: 10S 20M5I5D20M 10S

Is this a good way of defining the formula?

pid1 := 100 * M / (M+N+I)

or maybe

pid2 := 100* M / (M+N+I-D)

the example: pid1 = ~88% but then pid2 = 100% (which wouldn't make much sense).

Related: http://biostar.stackexchange.com/questions/9358/is-there-any-r-package-to-parse-cigar-element-of-sam-files/17031#17031

Thank you very much

Answer: No one seems to have mentioned it, but how about the "NM" tag ("edit distance" in SAM format specification)? I think this is the closest to "mismatches" that can be obtained in a relatively fast way by parsing bam files. One caveat, indels are always going to be a problem when defining identity and I don't know how they are counted in the "NM" tag, perhaps different programs have different ways of counting them.
Question: Hi, first post here.
So I'm trying take the CDS out of various species' orthologous sequences. I'm running on a Linux server, and am mainly aiming to use BioPython or Linux programs for this.

I've run OrthoFinder on 28 species of seaweed, which gave out roughly 10,000 orthogroup sequences fasta files, each of which is a a multi-fasta file. I've concatenated each of them into one huge multifasta file, and now I want to extract the fasta files according to their species into a new multifasta file (so 10k files -> 1 file -> 28 files, one per species).

How do I do this? I'm still fairly new to BioPython, so I'm still wrapping my head around things. I know I'll definitely need SeqIO, not sure what other libraries I'll need. I already have a text file with all the species listed, one per line.

Thanks heaps for any help.
Lachlan

Answer: Here is a Python code that reads one FASTA file and creates multiple FASTA files for each species separately.

from Bio import SeqIO

d = {}
fh = open('OG0000036.fa')
for seq_record in SeqIO.parse(fh, 'fasta'):
species_name = seq_record.id.split('-')[-1]
if species_name not in d:
d[species_name] = open(f"{species_name}.fa", 'w')
d[species_name].write(seq_record.format("fasta"))
fh.close()

**OUTPUT:**

FASTA files:

Ahnfeltiopsis_flabelliformis.fa
Betaphycus_philippinensis.fa
Ceramium_kondoi.fa
Chondrus_crispus.fa
Chroodactylon_ornatum.fa
Dumontia_simplex.fa
Eucheuma_denticulatum.fa
Glaucosphaera_vacuolata.fa
Gloiopeltis_furcata.fa
Gracilaria_blodgettii.fa
Gracilaria_lemaneiformis.fa
Gracilaria_sp..fa
Grateloupia_catenata.fa
Grateloupia_filicina.fa
Grateloupia_livida.fa
Grateloupia_turuturu.fa
Heterosiphonia_pulchra.fa
Kappaphycus_alvarezii.fa
Mazzaella_japonica.fa
Neosiphonia_japonica.fa
Porphyridium_cruentum.fa
Porphyridium_purpureum.fa
Pyropia_yezoensis.fa
Rhodochaete_parvula.fa

For example, `Grateloupia_filicina.fa` has three sequences:

>scaffold-ZJOJ-2006011-Grateloupia_filicina
MVVVFGLTGNDFALVVADMTSARSIMCFKHDEDKIERIDERKVLATAGEHSNRIEFSEYI
QKNLALMKLQTGLELSNHGTANFIRNEVAKALRTRGAYNTNSIMAGFDETGPAQKVNFTA
HGYASYFSLSVMDSKWRQDMTLEEGKKLVQECIDQLKSRFLINQPKFMMKIVTDQGITE
>scaffold-ZJOJ-2006903-Grateloupia_filicina
MDTLLGIAGEGFVVLAADAQVARSILLYKNDMDKIAHLSENKALACAGPQSDCVSFTEYI
SKNMALYELNNDVKLSTKAAASFIRGELAKALRKGPFQTQILMGGVDKRAAAEAEGKDDA
SLFWLDYLGTLQKVPYGAHGYGAAFTLSVMDREYVKGLSLDEALAIIDNCIKELHTRFLI
AQKNFVIKVVTAEGIK
>scaffold-ZJOJ-2055484-Grateloupia_filicina
MESIIALEGDGYVLIAADVASARSVVVMKDDMDKIRPLDSHKLFAAAGIPGDVSKFTEHV
QKDVRLYNMRSGITMSTAAAANYTRGELARFLRRSPFQCNVLIGGYDAPPYGHGPSLYSC
DYLGTLTKLKFAAEGYAQYFVLSTLDRYWKKNLSIEDGVEVIKKCVAEVQKRLVINQPRF
AIKIVDKNGVRVID
Question: Dear Friends,

My vcf file has SNPs available for different population(Africa, America, Europe,East Asia and South Asia ). I want to extract the data for Europe and East Asia together . Kindly let me know the possible ways.

Thanks in Advance

Answer: You can do this easily using vcftools, GATK tools, plinkseq etc.

you first have to generate a text file with the list of samples that form the population of your choice, let's say "population_of_interest.txt"
Then,

vcf-subset -e -c population_of_interest.txt input.vcf > output.vcf

or

vcftools --vcf input.vcf --keep population_of_interest.txt --recode > output.vcf
Question: Hi,

maybe it's a stupid question, but, are gene symbols supposed to be case sensitive or case insensitive?

E.g. GENE1, gene1, GenE1 represent the same gene?

Thanks

Answer: It'll depend on the organism. Human genes, for example, are always upper case (e.g., "GENE1"), with the exception of ORFs (e.g., C17orf12, which probably doesn't exist). Mouse genes, on the other hand, have a capitalized first letter and the rest lower case (e.g., "Gene1", which also needs to be italicized). In neither of those would a gene named "GenE1" exist. Having said that, you may see things like "GenE1" in print, since it's sometimes easier to see that that means gene family "Gen" subgroup "E" member "1", though it's often not technically correct. Technically a drosophila gene could be named "GenE1", since there's no strict rule on case (everything after the first letter is normally lower-case, but there's no "must" in that rule).
Question: Hi All,

Is it possible to convert a bam index (.bai) to human readable format? I'm asking mostly for learning, I don't plan to do anything special with it.

There is a good description of the bam index in the [sam spec][1] document, session 5.2 but I would like to see a "real" example of an index in plain text.

Thanks

Dario

[1]: https://samtools.github.io/hts-specs/SAMv1.pdf

Answer: You can write such a tool yourself using R (just as an example):

> x = file('/Users/sdavis2/Downloads/vcfanno_0.0.7_darwin_386/example/ex.bam.bai')
> open(x,'rb')
> readChar(x,4) # Magic number
[1] "BAI\001"
> readBin(x,integer(),size=1,signed=FALSE) # Number of Ref Sequences
[1] 86

At this point, you can use loops to loop over each of the blocks, reading in data as described in the spec.
Question: Hello,

this is my first post here. I would like to know if there is a way to call IGB (Integrated Genome Browser) through Java. I want to start a specific instance of IGB from a Java application and automatically show some tracks of interest that are on local file. Is this possible? Is there a way to use IGB into a Java application? or at least is there a way to use a whatever genome browser like that?

Thanks in advance

Answer: Hi,

Thank you for the great question!

Here are some additional possibilities that may be a good fit for what you want to do.

**Option One: Use IGB links (probably the simplest) - port 7085**

One option is to use ports on localhost, as IGB "listens" to port 7085 on localhost and can respond to commands encoded in a URL.

Here is some documentation about this to get you started:

- https://wiki.transvar.org/display/igbman/Controlling+IGB+using+IGB+Links
- https://wiki.transvar.org/display/igbman/Scripting+and+the+IGB+command+language

**Option Two: Use IGB commands - port 7084**

Starting with IGB 8.2 (release next week!) you can control IGB by sending commands directly to port 7084.

We implemented this to enable controlling IGB from within an interactive R session, but you could use this same mechanism to control IGB from other applications implemented in other languages.

To see how the R code looks, look at SRAdb/IGBFunctions.Rmd in https://bitbucket.org/lorainelab/rtoigb

Note that this is based on similar functionality in Bioconductor package SRAdb, hence the name.

**Option Three: OSGi services (requires Java coding)**

Possibly another option would be to use IGB's OSGi plug-in architecture to hook directly into IGB functions that are exposed as OSGi services. That would require more coding -- but if you wanted to invest the time, you could probably do some very nice stuff.

To find out more about how to work with IGB as a developer, see the developer's guide:

https://wiki.transvar.org/display/igbdevelopers/Home

Thank you again for your great question. I will post again here as we make more demonstration code available.

Yours,
Ann
Question: Im using this code to make based on log2foldchange and padj value ,im getting the plot but i want those value for my reference how do i extract the same .

alpha <- 0.05 # Threshold on the adjusted p-value
cols <- densCols(res$log2FoldChange, -log10(res$pvalue))
plot(res$log2FoldChange, -log10(res$padj), col=cols, panel.first=grid(),
main="Volcano plot", xlab="Effect size: log2(fold-change)", ylab="-log10(adjusted p-value)",
pch=20, cex=0.6)
abline(v=0)
abline(v=c(-1,1), col="brown")
abline(h=-log10(alpha), col="brown")

gn.selected <- abs(res$log2FoldChange) > 2.5 & res$padj < alpha
text(res$log2FoldChange[gn.selected],
-log10(res$padj)[gn.selected],
lab=rownames(res)[gn.selected ], cex=0.4)

when i view gn.selected i get only logical value that is true or false

Any help or suggestion would be highly appreciated

Update I'm doing this

> DF <- DF[DF$log2FoldChange > 1.5 & DF$padj < 0.05,]

is that suffice and am i doing it correctly ?

Answer: **Edit (October 24, 2018):**

This is now a Bioconductor package: https://www.biostars.org/p/335751/

---------------------------------------------
---------------------------

Your code appears to run fine on my DESeq2 results objects:

-----------------------------------

I normally do these simple volcano plots a different way:

par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)

topT <- as.data.frame(resultsObject)

#Adjusted P values (FDR Q values)
with(topT, plot(log2FoldChange, -log10(padj), pch=20, main="Volcano plot", cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))

with(subset(topT, padj<0.05 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(padj), pch=20, col="red", cex=0.5))

#with(subset(topT, padj<0.05 & abs(log2FoldChange)>2), text(log2FoldChange, -log10(padj), labels=subset(rownames(topT), topT$padj<0.05 & abs(topT$log2FoldChange)>2), cex=0.8, pos=3))

#Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05
abline(v=0, col="black", lty=3, lwd=1.0)
abline(v=-2, col="black", lty=4, lwd=2.0)
abline(v=2, col="black", lty=4, lwd=2.0)
abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col="black", lty=4, lwd=2.0)

--------------

There is an even better solution that I and colleagues developed using *ggplot2*, which allows you to easily fit labels into your plot using `ggrepel()`.
Question: bgzip files are backward compatible with gzip, but I have issues when using bgzip compressed vcf files with snpeff (java) or perl scripts that uses IO::Uncompress::Gunzip (that I believe it uses zlib under the hood). In both cases the data is decompressed but truncated after a few hundred lines aprox. I could be totally wrong but I was wondering if zlib (or whatever gzip compatible library they are using) is getting confused with the bgzip bloks and only processing one or a few of them leaving the output incomplete.

perl code that does not work:

#!/usr/bin/env perl
use strict;
use IO::Uncompress::Gunzip qw(gunzip $GunzipError) ;

my $infile = shift;
my $infh = IO::Uncompress::Gunzip->new( $infile ) or die "IO::Uncompress::Gunzip failed: $GunzipError\n";
my $line_count = 0;
while (my $line=<$infh>){

$line_count++
}
print "total lines read = $line_count\n";

This gives 419 lines

$ perl /home/pmg/tmp/test_zlib-bgzip.pl 460112_TTAGGC_L005_L006_C3HVJACXX.sorted.rmdup.varsit.vcf.gz
total lines read = 419

but using open with gzip pipe works:

#!/usr/bin/env perl
use strict;

my $infile = shift;
open(my $infh , 'gzip -dc '.$infile.' |'); # I can use bgzip intead gzip
my $line_count = 0;
while (my $line=<$infh>){

$line_count++
}
print "total lines read = $line_count\n";

Gives the expected number of lines

$ perl /home/pmg/tmp/test_gzip-bgzip.pl 460112_TTAGGC_L005_L006_C3HVJACXX.sorted.rmdup.varsit.vcf.gz
total lines read = 652829

I googled about and I was unable to find quickly any relevant entry, but this is something that I am sure other people would have already faced. Do someone have a clue about why is this happening? I am using ubuntu 12.04.4 with perl 5.16

Answer: Hi. I have just encountered the same problem and it seems there is a pure-Perl solution surrounding this as shown below

my $status = gunzip("in.gz" => "out.gz", MultiStream => 1)
or die "gunzip failed: $GunzipError\n";

It also seems Apache Commons Compress can provide similar support for Java. Hope this helps
Question: Hi all,

Sorry I know that this question has been asked several times, but unfortunately I haven't been able to find the right answer, or didn't understand.

I'm trying to get TMM normalized counts thanks to edgeR.

I understand that I have to compute normalization factors :

dgList <- calcNormFactors(dgList, method="TMM")

which gives me a normalization factor for all samples :

head(dgList$samples)

group lib.size norm.factors
S1 1 21087314 0.9654794
S2 1 16542810 1.1589117
S3 1 18875473 0.8763291
S4 1 15865414 1.0864038
S5 1 19179795 1.0488230
S6 1 15063992 1.0707007

But at this step I don't know what to do to get a matrix of normalized TMM counts.

I know that I can get CPM normalized counts thanks to :

cpm(dgList)

But CPM and TMM are not the same, right ?

Thanks in advance for any of your input on this topic.

Answer: No, CPM and TMM are not exactly the same indeed.

perhaps try this snippet of code:

dgList <- estimateCommonDisp(dgList)
dgList <- estimateTagwiseDisp(dgList)
norm_counts.table <- t(t(dgList$pseudo.counts)*(dgList$samples$norm.factors))
write.table(norm_counts.table, file="./normalizedCounts.txt", sep="\t", quote=F)
Question: We are currently downloading and analyzing multiple large WGS datasets (30-50x) from patients. So far, we downloaded the data of 20 patients from dbGaP/NCBI (tumor and matched normal respectively). More samples are planned to be included. The download itself via prefetch/fasp was relatively fast and smooth but now the problems begin, so maybe you have some experience in how to optimize things.

-Sra to fastq via fastq-dump is often unbearably slow. Not only is the fastq-dump slow itself, but I often experience I/O bottlenecks on our university cluster, which uses gpfs (not lustre as I stated yesterday). fastq-dump is often stuck in "D"-state, so uninterruptable sleep.
To speed up things, I dumped large sra into several fastq files, using the -N and -X options, but merging these chunks via GNU cat was also extremely slow, sometimes with only a few hundred MB in several hours. Is that normal (the server does not run on SSD as far as I know).

-same goes for alignment sorting. I tried to rather use fewer threads with SAMtools sort but more memory per thread to avoid the creation of too many tmp files, which then need to be merged again. still, even merging few (< 50) files takes hours and hours of time, again with only few hundreds MB in several hours. That often collides with the walltimes.

It would be great if you could share your experiences in how to handle these Terabyte-scale data, and what tricks one can apply in order to avoid performance bottlenecks.

UPDATE: It seems that the main bottleneck is reading the files from disk, rather than writing them after being processed.

Answer: The solution we came up was the following: Our file system is simply slow, and there was nothing that could really be done about it. The main bottleneck was reading from the file system, rahter than writing. Fortunately, some of the nodes had local SSDs, which I could use. So loaded the SRAs via prefetch (ascp) to the SSD, then fastq-dump them from there, outputting directly to /scratch. Thanks to ascp, the download of a 40-100Gb files was done in no time, and the dumping was speeded-up by (never benchmarked it) I think factor 10. Thanks very much for all your suggestions.
Question: Dear all,

I have indexed bam file, and I want to print all alignments on chromosome "2" using pysam.

My code is

import pysam
bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
for line in bam.fetch("2"):
print line

What I get is:

D00733:162:CADM2ANXX:2:1311:7466:63232	163	18	3328	255	1S47M2S	18	3419	47	TTGCTTAGTGTCCGAAATACCATCCTCAAGGCTAAGAACTAAATCGATTA	array('B', [33, 33, 33, 33, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]
D00733:162:CADM2ANXX:2:1215:11810:89291	163	18	3356	255	44M6S	18	3441	44	GGCTAAGAACTAAATCGATTATTCTGGCTCGTAACGCATATAATATGGCA	array('B', [33, 33, 33, 32, 33, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 35, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 37, 38, 38, 38, 38, 38, 38, 33])	[('NH', 1), ('HI', 1), ('AS', 70), ('nM', 10)]

So there are 2 questions:
1. Why pysam does not preserve format of bam file?
2. Why does it print out chromosome 18 instead of 2?

I would appreciate if somebody points out on what am I doing wrong.
Thanks

Answer: Hi, you are not wrong, you are doing it correctly :)

The AlignedSegment object is independent of a SAM-file. The AlignedSegment contains only an index of the chromosomal identifier (tid), and, thus, it does not contain the chromosome number information. It means that the number you see as 18 in the field 3 of the line is the number corresponding to chromosome 2 in the sequence dictionary. You can check that by doing the following:

import pysam
bam = pysam.AlignmentFile("Aligned.sortedByCoord_rep1.out.bam", "rb")
for line in bam.fetch("2"):
print line
print line.reference_id
print bam.get_reference_name(line.reference_id)

As you might have notice, I have added two lines to your previous code:

1.`line.reference_id` contains the index of the reference sequence in the sequence dictionary

2.`bam.get_reference_name(line.reference_id)` will return the reference name corresponding to numerical tid. In your case chromosome 2.
Question: Hello guys
I'm working on Rna SEQ data, and I want to get good conceptual knowledge about the analysis of the process.

I used the tool bowtie for indexing my reference genome and then the rest is according to the TUXEDO pipeline.

I don't understand what exactly, "indexing" of the reference genome is and why it had to be done.

I do get the point that there are many genes that can be identified and can be matched to the genome and indexed, But i am not clear as ti why should we do this , how is it done . ( i am new to bioinformatics and with minimal knowledge its so hard to understand the algorithm part of it )

so can anyone put it in simple concepts .

Answer: Indexing a genome can be explained similar to indexing a book. If you want to know on which page a certain word appears or a chapter begins, it is much more efficient/faster to look it up in a pre-built index than going through every page of the book until you found it. Same goes for alignments. Indices allow the aligner to narrow down the potential origin of a query sequence within the genome, saving both time and memory.

You may have a look at Ben Langmeads youtube channel (creator of the bowtie and bowtie2 alignment software), where he lectures about concepts of sequence alignments, including the idea of indexing.

https://www.youtube.com/user/BenLangmead