Standard for when to collapse phylogenetically uninformative nodes I have been creating phylogenetic trees using RAxML (using Boostrap + ML) and I end up with a bestTree file, outputted in newick format. When I open up the tree in FigTree, I notice that it collapses many of the internal nodes with small values (1 x 10-6). I am working on visualizing the trees using ETE, which when I render the tree includes these small breakpoints. I was wondering if there is a standard for when to collapse nodes that are considered phylogenetically uninformative that I can apply. Thank you for any help.

If 1E6 is genetic distance, this means one nucleotide difference per million base pairs, which is of course incredibly small and bootstrapping at this level doesn't really work. So what I suspect is happening is that the bootstrap is so low it is not present on the tree. You can assign a cut-off bootstrap but anything below 75% is not interpretable in any case. What might be happening is you are producing the bootstrap "cladogram" which is showing a default branch length, but when <1E6 genetic distance it is showing a trifurcation. This simply means the bootstrap of the consensus tree is so low it is collapsing the branch and this is a function of the very small genetic distance at this node.

Looking at your treefile I don't see any collapsed clades.
(Group_9:0.00000100000050002909,(Group_50:0.00507659947127924935, (.

The clade above has not been colapsed it is just a very small distance in a paraphyletic group. It looks fine to me, i.e. Group_9 is 1E6. It could look 'collapsed' simply because it is so small it is difficult to see. I hope thats ok?


DNA bioinformatic search and methylation analysis by qMSP Im quoting the methods which was used

  We performed a bioinformatics analysis to identify potential CpG
  islands in the genomic sequences of RANKL and OPG genes. The USCS
  Human Genome Browser public database (genome.ucsc.edu), the Primer
  Express software (available from Applied Biosystems), and CpG Island
  Explorer software (cpgie.sourceforge.net) were used to explore the
  nucleo- tide sequences surrounding de transcription start site (TSS)
  of RANKL isoforms I and II (NG_008990.1 Ref Seq Gene), and OPG
  (NG_012202.1 Ref Seq Gene). A (C + G)/total bases ratio

    0.5 and CpG observed/CpG expected ratio >0.6 were used as criteria.

Im looking for the region which they have used

This is the paper

"Role of DNA methylation in the regulation of the
RANKL-OPG system in human bone"

I tried using ucsc to check the region but didn't get much. The idea is to find the region to design the primer for the same

How can i get the region from ucsc ,the CpG island around the TSS of the above mentioned genes which are RANKL and OPG

Any suggestion or help would be highly appreciated

It's likely that they accomplished this by using the UCSC Table Browser to:

Find the loci of the two isoforms you mention
Finding CpG islands (Group = Regulation, Track = CpG Islands) near those transcription start sites
Design primers using Applied Biosystem's software

Alternatively, they list the exact primers used in Table 1.

They don't detail exactly how they did it, so if you have specific questions it might be worth it to email the corresponding author


How to incorportate RIN values as covariate in the design matrix? I have been following the last DESeq2 pipeline to perform an RNAseq analysis with a dataset with low rin samples in the experimental (or treated) and high rin on the control ones.

I read a paper in which they perform RNAseq analysis with time-course RNA degradation and conclude that including RIN value as a covariate can mitigate some of the effects of low RIN in samples.

My question is how I should construct the design in the DESeq2 object:

~conditions+rin
~conditions*rin
~conditions:rin

none of them... :)

I cannot find proper resources where explain how to construct these models (I am new to the field...) and I recognise I crashed against a wall with these kinds of things. I would appreciate also some links to good resources to be able to understand which one is correct and why.

Thank you in advance!

In the DESeq2 manual there is a section titled "How can I include a continuous covariate in the design formula?" that deals with your question.
Basically the process is no different from using a discrete covariate.

DESeq <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ RIN + condition)

In this example condition is the experimental groups and RIN is the continuous covariate for the RNA-integrity.

For this to work properly there has to be a linear relationship between the RNA integrity and the gene expression values. As stated in the manual: "Continuous covariates might make sense in certain experiments, where a constant fold change might be expected for each unit of the covariate."
Therefore, it may be better to code the RIN variable as a factor like RIN_high or RIN_low. This is more robust than using the absolute values of the variable itself.

However, from your desription I get the impession that the RIN values are generally lower in the treatment group compared to the control group. If this is the case, using RIN as a covariate will by of limited use. DESeq2 will not be able to discriminate if a expression difference is due to a difference in RNA integrity or due to your treatment. Also, in this case coding RIN into high and low would result in an error because the RIN and treatment variable would define the same groups.

Generally, covariates work best when they are evenly distributed across the experimental groups.


Assemble reads for a specific gene I have a lot of unassembled sequencing data and I want to build phylogeny for some genes from these data. I can retrieved the sequene by mapping reads to a reference sequence but if the two species are too divergent, reads will not match and I won't get the full CDS. Is there a better way to deal with that? Maybe I can do that to retrieve the well conserved part and then extend this region with reads that overlap (as if it was gap-filling)?
What do you think of this approach? Do you know a software for such extension process?
Thanks

That's a good question. I remember seeing a poster about making trees directly from reads, but unfortunately, I can not find any other record than this tweet. Perhaps you can try to email David, if he would be keen to share the implementation.

While searching for any information about David's work, I run into a GitHub repository called reads2trees. It's quite old and it's not super popular, so hard to say if it can be useful. But the name matches the problem description quite well.

I also remember there were some attempts (An example paper) to get phylogeny directly from kmers, but I have never tried that on my own. I also think these approached were focused on microbial phylogienies. So again, not sure if it's going to be useful for you.

That's all I got.


SNP located within a promoter region (pig) I have a couple of SNP identifiers such as MARC0073381 or ALGA0066960. The corresponding platform is Illumina Porcine SNP60 BeadChip (WG-410).

I want to know if these SNP are located within a promoter region of a gene.

Bioconductor offers information for Affymetrix chips, but I am not able to map the identifiers.

How can I map SNPs from the beadchip to the porcine genome?

As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if there are only a few locations of interest it would be easier to use grep than try to work out the right way to load it into R (or similar).

For completeness, you want the lines between [Assay] and [Controls]. Here's an example way to pre-process the "CSV" file into a more standard format that just has one data table in it:

$ grep -n '^\[' PorcineSNP60v2_15031945_C1.csv
2:[Heading]
7:[Assay]
61574:[Controls]
### See that the table starts at line 8 (after [Assay]),
###   and finishes at line 61753
$ tail -n +8 PorcineSNP60v2_15031945_C1.csv | head -n $((61574-8)) > cleaned_PorcineSNP60v2_15031945_C1.csv

But for a few markers and grep, it doesn't matter if you're working with the clean or dirty CSV file:

$ grep -e '^IlmnID' -e '^MARC0073381' -e '^ALGA0066960' PorcineSNP60v2_15031945_C1.csv
IlmnID,Name,IlmnStrand,SNP,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,AlleleB_ProbeSeq,GenomeBuild,Chr,MapInfo,Ploidy,Species,Source,SourceVersion,SourceStrand,SourceSeq,TopGenomicSeq,BeadSetID,Exp_Clusters,Intensity_Only
ALGA0066960-2_B_F_2199219285,ALGA0066960,BOT,[T/G],0054749354,TGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA,,,10.2,12,55218591,diploid,Sus scrofa,PorcineSNP60v2,2,BOT,GGCTCACCTCTGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA[T/G]CCAGGCTTCTCTCTCCTACTCCTAGGGCCCCTGATGGTTCCTGCATCCTGACCAATAGTG,CACTATTGGTCAGGATGCAGGAACCATCAGGGGCCCTAGGAGTAGGAGAGAGAAGCCTGG[A/C]TAGACCTGTATCTTGCGAGTGTTGGGCTTGAGCTGAGCAGACCTGTGGCAGAGGTGAGCC,670,3,0
MARC0073381-2_T_F_2199252059,MARC0073381,TOP,[A/G],0047761415,TGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA,,,10.2,7,61808556,diploid,Sus scrofa,PorcineSNP60v2,2,TOP,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,660,3,0

I can see from this that the first marker is Chromosome 12, position 55218591, and the second marker is Chromosome 7, position 61808556. The Illumina annotation doesn't include the nearest gene, so you'll need to hunt for those locations on a genome browser (or using mapping with Bioconductor).


Using sed:

sed 's,\([01][/|][01]\):[^ ]*,\1,g' file
#    ^  ^               ^     ^  ^
#    |  |               |     |  *--- applies to all found pattern
#    |  |               |     *------ print the content inside the parenthesis
#    |  |               *------------ any character except a space
#    |  *---------------------------- the pattern to be kept
#    *------------------------------- substitution

You also may want to format the output using column command:

sed 's,\([01][/|][01]\):[^ ]*,\1,g' file | column -t
#CHROM  POS     ID  REF   ALT   FATHER  MOTHER  DAUGHTER1  DAUGHTER2  DAUGHTER3   SON1  SON2
1       736689  .   T     A     0/0     1/0     0|0        0|1        0|1        0|0   0|1
1       758663  .   CGCC  TTAG  0/0     0/1     0|0        0|1        0|1        0|0   0|1

Note the comma separators in the sed command can replaced with other characters typically "/"


Effect of mutation in DNA sequence on transcription factor binding sites How much does a single mutation/alteration of a nucleotide affect the presence of a transcription factor binding site (TFBS)?
I am from computer science background(Obviously).
I want to make a general assumption about the number of mutated bases in a TFBS to minimize the binding affinity?

I haven't worked on TFBSs but I really, really doubt you will be able to do this. In mutations, like in real estate, it's location, location, location. It isn't a question of how many mutations but what mutations and where. TFBSs are usually described using position weight matrices, precisely because that's how you can show the importance of different changes in different positions.

Here are a couple of random Google image search results for "transcription factor binding sites logo":

(This was taken from here and you might want to read through it.)

(source: Front. Genet., 23 February 2016 | https://doi.org/10.3389/fgene.2016.00024)

The main take home messages from these images are:

Different positions have different importance and potential to disrupt binding. In the first image, I would expect that changing the T at position -3 to anything else would be pretty bad. However, changing an A to a T at -8 would likely be neutral.
Different TFs have very different binding sites and very different binding site lengths.  Look at panel B of the second image. That one appears to be a very conservative site and likely any change would be deleterious. Conversely, changing any of the nucleotides of the 17bp spacer shown in panel A would likely not be important as long as you don't change the length.

Please note that I have not done more than glance at the sources I took these from, so my inferences may well be wrong, I only use them to illustrate the general point: you cannot really put a number to this. You will always need to look at each specific change and each specific position. In some cases, a single change can disrupt binding while in others multiple changes might have no effect at all.

You might want to look into TFBS databases like JASPAR or TRANSFAC.


Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset? As far as I understood, for most assembly programs, the scaffolding step takes into consideration paired-end information in order to get from contigs (contiguous sequences) to scaffolds (longer sequences that might have some N-filled gaps).

My assembly software of choice, MEGAHIT, uses paired-end information to build the contigs, but it does not output a standard scaffold. So I am wondering the following, related things:

-is it meaningful to run a scaffold program on the output of MEGAHIT? I imagine there might be some instances in which paired-end information could span a gap.

-Which software would you recommend for it? (I've tried soapdenovo2 and SSPACE, but they appear not to be actively maintained so I have the issue of 'it doesn't work and I can't do anything about it')

-Could relevant information (eg. two contigs being connected by pair-end information) be recovered by alternative and perhaps more user-friendly means, such as exploring the assembly graph with Bandage?

Thank you for your time!

Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.

I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.

However, I think that it is perfectly ok to scaffold using orthogonal data. Here are some examples of orthogonal data:

A reference assembly (weird for metagenomics but still possibly ok).
long reads, e.g. ONT or PacBio (low coverage probably fine).
Proximity ligation data, e.g. Hi-C reads.
Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
Linked-reads (e.g. 10X).
Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.

Full disclosure: I work for a company that sells Hi-C kits for such applications.

Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.


Get & Annotation of isomiR (isoform of miRNA) profile from RNA-seq I should get and analyse isomiR (isoform of miRNA) profile from RNA-seq with FASTA format. The best tool for this is "DeAnnIso". But here is a problem with this tool, it only allows you to upload FASTA format files  with a maximum size of 50MB. However, the size of my FASTA file is more than 1000MB. I tried to compress the file in different ways and uploading, but it doesn't work.
Does anybody know any way to solving this problem?

Or Does anybody know any other easy way to get isomiR profile from RNA-seq with FASTA format?

A sample of RNA-seq like this:

>RR42357.95441.1 I-D00105:274:H7MDXX:1:112:8124:865 length=51
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACCTTGTAATCTCGTATGCC
>RR42357.95442.1 I-D00105:274:H7MDXX:1:112:8133:866 length=51
AAAGTGCTACTACTTTTGAGTCTTGGAATTCTCGGGTGCCAAGGAACTCCA
>RR42357.95443.1 I-D00105:274:H7MDXX:1:112:8144:8637 length=51

Thank you very very much!

The obvious answer is to break up the file size and upload each separately,
First I'd check the file size,

wc infile.fa

Output

  XXXXX      YYYY     ZZZZZ infile.fa

The number you want is XXXX - the number of lines in the file.

You'd then decide how many files you wanted to break it up into ... in your case divide XXXXX by 3. Then .. for example

sed -n '1,5000000p' infile.fa > outfile.fa
sed -n '5000001,2000000p' infile.fa > outfile2.fa
sed -n '2000001,4000000p' infile.fa > outfile3.fa

I'm not certain this would be cool however. An miRNA expert could answer this question very quickly. I think is you should check your data for sequence duplication. Simply piping the fasta file through a simple distance matrix and screening for "0" values would confirm this.

Creating a frequency table of duplicate RNA must be a feature of Bioconductor, but in any case screening for duplication (technically duplicate records) is a formal SQL database command.

It can be done a Perl/Python but its a bit clunky for a data set this size. I'd do it by setting the sequence as the key within a hash, the inefficiency has been solved (but its a weird way to do it).

You would have a database for "headers" for a given sequence, which would take up space, but I would assess the potential level of duplication first.


How to specifically select cluster on Featurescatter on Seurat? Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger.

FYI, The filtered_feature_bc_matrix includes an aggregated population of tumor-infiltrating immune cells that are derived from 6 different syngeneic mouse tumors (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861). I downloaded fastq files and ran "cellranger count" and "cellranger aggr" commands to obtain combined data.

My strategy includes the use of the "feature scatter" function deposited in the library(Seurat). Finally, I ended up getting results below.

My question is that can we select a specific cluster(s) or identity on this plot? For example, I want to see correlation values in cluster 2 only.
The usage of Feature scatter is as follows.

FeatureScatter(
  object,
  feature1,
  feature2,
  cells = NULL,
  group.by = NULL,
  cols = NULL,
  pt.size = 1,
  shape.by = NULL,
  span = NULL,
  smooth = FALSE,
  combine = TRUE,
  slot = "data"
)

Do you know how to get the Pearson correlation of my two genes using Seurat? This is the best I can do.

Thank you.

You can use the cells argument to only display the data for a specific subset of cells.

For example to extract the values for a cluster 2:

# The cluster identity is stored in the meta.data
# In this example there is a clustering column called "RNA_snn_res.0.7"
# but any metadata column can be used

clust2.cells <- rownames(Seurat.data@meta.data[Seurat.data$RNA_snn_res.0.7 == 2, ])

FeatureScatter(Seurat.data, feature1 = "gene_A", feature2 = "gene_B", cells = clust2.cells)


What is the --rev_comp_mapping_barcode parameter in the QIIME1 script extract_barcodes.py? I'm using the QIIME1 scripts extract_barcodes.py to extract barcodes from a dual-indexed MiSeq library, and then demultiplex my reads using split_libraries_fastq.py.  I have been having trouble correctly using a couple of parameters relating to reverse complementing barcodes.

When extracting the barcodes with extract_barcodes.py, I use --rev_comp_bc1 (as directed by the individual in charge of sequencing). Then when I demultiplex the reads with split_libraries_fastq.py, sometimes I need to include the parameter --rev_comp_mapping_barcodes to get the script to work, and sometimes I don't.  I'm not sure why, unless the barcodes I'm receiving in my mapping files are sometimes written in different orientations (I don't write them).

I think the parameter --rev_comp_bc1 is self-explanatory, but could someone explain what --rev_comp_mapping_barcodes does and when to use it?
To make things simple we can assume that the two barcodes are written from 5' to 3' in the mapping file.

As the help page says:

  --rev_comp_mapping_barcodes
      Reverse complement barcode in mapping before lookup (useful if barcodes in mapping file are reverse complements of golay codes) [default: False]

As the help page of extract_barcode says:

  --rev_comp_bc1
  Reverse complement barcode 1 before writing [default: False]

You seem to be doing two complete different things. extract_barcode.py is designed to

  format fastq sequence and barcode data so they are compatible with split_libraries_fastq.py

While the first one is for:

  demultiplexing of Fastq sequence data where barcodes and sequences are contained in two separate fastq files

So it seems that depending on the machine used or the previous processing you are given with the barcode in another file or together. I would advice to ask the person who gave you this data.


How to subset a VCF by chromosome and keep the header? I would like to subset a VCF which only has chromosome 2.

The problem with using various grep commands, e.g.

grep -w '^#\|^2' my.vcf > my_new.vcf

or if there's a 'chr' prefix

grep -w '^#\|^chr2' my.vcf > my_new.vcf

is that this will remove the header.

I've been trying to use vcftools with the following command:

vcftools --vcf input.vcf  --chr 2 --out subset

which does not output a VCF subset.vcf as expected, but rather subset.log:

VCFtools - 0.1.15
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf input.vcf
    --chr 2
    --out subset

After filtering, kept 2 out of 2 Individuals
After filtering, kept 80 out of a possible 720 Sites
Run Time = 0.00 seconds

Given that the run time is Run Time = 0.00 seconds, I guess there's an error here.

How can I subset a VCF and keep the header? Is there an option with bcftools instead maybe?

This seems to work for me with bcftools filter and the -r or -R argument.

   -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]

  Comma-separated list of regions, see also -R, --regions-file. Note that -r cannot be used in combination with -R.

   -R, --regions-file FILE

  Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the
             tab-delimited file are:
             CHROM, POS, and, optionally, POS_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while
  bgzip-compressed and
             tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note
  that chromosome
             ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE.
  However, within chromosomes, the
             VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that
  overlapping regions in FILE can
             result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be
  used in combination with -r.

Here's an example:

$ bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20171222
##source=PLINKv1.90
##contig=<ID=1,length=249167691>
##contig=<ID=2,length=242695901>
##contig=<ID=3,length=197800245>
##contig=<ID=4,length=190915651>
##contig=<ID=5,length=180666277>
##contig=<ID=6,length=170877445>
##contig=<ID=7,length=159086805>
##contig=<ID=8,length=146293415>
##contig=<ID=9,length=141018424>
##contig=<ID=10,length=135434552>
##contig=<ID=11,length=134938471>
##contig=<ID=12,length=133763353>
##contig=<ID=13,length=115045730>
##contig=<ID=14,length=107285438>
##contig=<ID=15,length=102369712>
##contig=<ID=16,length=90141356>
##contig=<ID=17,length=81006630>
##contig=<ID=18,length=78014583>
##contig=<ID=19,length=59071322>
##contig=<ID=20,length=62906515>
##contig=<ID=21,length=48077813>
##contig=<ID=22,length=51156934>
##contig=<ID=23,length=154847490>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746
##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  1_125   3_5555  4_7565  5_9
4   71566   rs13125929  T   C   .   PASS    PR  GT  0/1 0/1 0/0
4   87113   rs2006748   T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   110646  rs11727494  C   T   .   PASS    PR  GT  0/1 0/0 1/1
4   142550  rs11735742  T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   200631  rs6826124   A   G   .   PASS    PR  GT  0/1 1/1 0/1
4   221623  rs7695945   G   T   .   PASS    PR  GT  0/0 0/0 0/0


Significant gene set testing - limma While I have used my title, it is more of a general statistics question regarding how these significance tests are carried out, these include decideTests and mroast / roast among others.

To give a brief overview of my experimental design, I am analysing bulk RNASeq datasets and there are two main factors in the design matrix: patientID and casetimeViremia; where the latter identifies whether a patient is viremia positive or negative at a particular time point (1-4). This leaves us with 5 different possibilities for casetimeViremia: 1-min, 2-plus, 3-min, 3-plus, 4-min.

The design matrix was created in R by running:

model.matrix(~patientID + casetimeViremia, data = dgeObject$samples)

Where dgeObject is an object of the function DGEList.

For the purpose of answering this biological question, we wish to conduct the differential expression testing by seeing whether all the previously mentioned casetimeViremia timepoints are separately and significantly expressed relative to each person's first timepoint, which is labelled as: 1-min.

Now, I understood that should I wish to conduct this test, I do not need to specify a contrast matrix in the mroast function and can just do this by calling the column number of interest.. But in that case, how does mroast know that I wish to do the contrast relative to the 1-min timepoint? And, just wondering, should I wish to define a contrast matrix, how would I do this given that there are two colnames missing from the design, in this case one patientID and the 1-min timepoint?

Thank you all for your time.

The mroast function has an argument to specify which contrast do you want to test, quoting from the help page:

  contrast
  contrast for which the test is required. Can be an integer specifying a column of design, or the name of a column of design, or a numeric contrast vector of length equal to the number of columns of design.

So if you have the contrast of 1-min timepoint positive vs 1-min timepoint negative you can use this.
Read the help pages in detail, specially the help pages in limma are quite extensive and detailed to answer most questions if not all about how to use the software.


How to extract gene expression tables from this GEO dataset? I've downloaded this GSE43013 dataset using GEOquery in R. My understanding is that it contains expression data from liver, kidney, and brain for multiple species. I would like to produce gene expression tables for each species and tissue type. The goal would then be to perform PCA and regression analyses to see if there is a relationship between expression of certain genes and e.g. taxonomic position, diet, etc.

However, I have never worked with gene expression data before so I have very little understanding of the GSE file. I have failed to extract any kind of expression data from the file. Solutions that I have found online do not seem to be working for this dataset (for example producing empty tables where I assume there would normally be data on expression levels). Because of my inexperience I cannot determine why they are not working in this situation.

I would appreciate some help in understanding the format of this file, but most importantly, in figuring out how to pull expression data from it, if that is possible.

EDIT: I still haven't managed to pull data from the GSE file, but looking through the samples associated with the data set (such as this one), I can see that there are files formatted like this:

mus.lv.3

ENSMUSG00000000001  26.3376

ENSMUSG00000000003  0

ENSMUSG00000000028  0.501533

ENSMUSG00000000031  0.116078

I assume that the values in the left column refer to specific genes. Perhaps the values in the right column represent expression intensity? If so, would these values be normalized and comparable between species as they are?

Assuming that you have downloaded processed file, according to your link, it is either plain text (compressed) or GTF. If you would have plain text you probably would have posted a snapshot or sample of it, so I assume you have the GTF. According to format description for GTF-files:

Format consists of one line per feature, each containing 9 columns of data, plus optional track definition lines.
The 9 columns (or fields) are tab-separated and all fields contain a value (missing/NA-values are usually denoted with ".").
Fields are:

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

You can find more information on it here and here.

IF you don´t have the file and just have the query output, you probably will find the answer here, where extraction of expression data has been discussed in more detail.


Single sample in group: normal pipeline or Kal's Z test As stated, which one is better for differential expression analysis? When I say normal pipeline I mean limma-voom, edgeR and DESeq2 pipeline for standard analysis. Kal's z test is mentioned in this paper and I know some pre-made tools like the one on QIAgen uses this as well for single vs single comparison

The situation I have is that I have RNA-Seq data for some patients with rare translocations in leukaemia. There are two of these translocation with almost 10 or even more patients so I have no problem with them. And then there are these translocation which mostly I have only 1, some with 2 or 3 patients, for 7 different translocations.

Our hypothesis is that the two larger group can certainly form their own subtypes, and the others may fall into one of the three cases: (1) belongs to one of the two groups despite of the difference in translocation; (2) belongs to a separate group on their own i.e. complete outlier; (3) belongs to a separate group with other members in the "others" i.e. forming a third subtype. And this can be tested with clustering.

For that, I need to do supervised feature selection to filter out useless genes. I can do that with those with multiple samples, but not those isolated cases. So I thought I can get some DGE analysis and use a certain adj.p and fold change as cut off, or adj.p as ranking matrix, and select those genes for clustering.

So my main questions are:
(1) when using limma/edgeR/DESeq2, should the model design include all the translocation labels in one go or should I just add the one contrast I want to make and build a new model everytime?
(2) Is Kal's z test useful for one vs all DGE analysis and is it better?
(3) Are there other outlier analysis working better under this situation?

Many thanks!

You'll have a bit more power if you use all of the samples at once.
Nothing is useful for this, don't waste your time with it.
You're not looking for outliers, you're looking to see how samples cluster. You can just use correlation for that.


What is "unmapped read segments" in the output of samtools idxstats? samtools idxstats produces a four column output  (see here)
The output is TAB-delimited with each line consisting of reference sequence name, sequence length,
# mapped read-segments and # unmapped read-segments. It is written to stdout. Note this may count reads
multiple times if they are mapped more than once or in multiple fragments.

What is the meaning of unmapped reads here? Naively, I would take it to mean all the reads that didn't map to a specific feature... however in a sample with millions of reads most of them should be counted as such, which is not the case in my calculations, where the number of unmapped reads is of the same order as the number of mapped reads.
I will very much appreciate explanations from people experience with samtools.

An unmapped read segment is any entry with bit 4 of the flag field set (i.e., it's unmapped). Assuming the reference to which you aligned reasonably matches the sample from which the reads were generated, less than ~10% of to total reads should be unmapped.
Please note that BAM files have no idea what a "feature" is. They only contain information about reads and the chromosomes/contigs to which they align.


Why is there a chloride ion in this 3D model? Here we have a 3D model of a sequence that match with this sequence
VVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIK
in Blast

the 3D structure is here

Why is there a chloride ion in this 3D model?

It is a crystallisation artefact. Namely, the protein are placed in a condition where they fall out of solution without aggregating a gloop (as happens in most well in xstal trial) and to pack as a crystal they need to be placed nice and orderly. It is coordinated by water so is not relevant. A lot of structures have these —DMSO and ions are the most commonly.

Two footnotes.

1ZV7 is a crystal structure, i.e. solved by X-ray diffraction. A model is generally used to described an in silico model, such as those from Swissmodel, phyre, I-Tasser etc. These are generated ab initio (from scratch) or threaded (mapped on top of) against crystal/EM/NMR structures. For Covid I-Tasser has a nice picture of its models.
The best source of info about structures is the PDB. RCSB PDB is more commonly used, while PDBe is more tricker to navigate (or used to be). Both have the same data, but are better for the extras than NCBI.


Normalization methods with RNA-Seq ERCC spike in? ERCC spike-in is a set of synthetic controls developed for RNA-Seq. I'm interested in using it to normalize my RNA-Seq samples. In particular, I'd like to use the spike-ins to remove technical bias and any variation that should not be part of my analysis.

The site doesn't give any details on how I can do that.

Q: What are the possible normalization strategies? Can you briefly describe them?

You may consider using RUVSeq. Here is an excerpt from the 2013 Nature Biotechnology publication:

  We evaluate the performance of the External RNA Control Consortium (ERCC) spike-in controls and investigate the possibility of using them directly for normalization. We show that the spike-ins are not reliable enough to be used in standard global-scaling or regression-based normalization procedures. We propose a normalization strategy, called remove unwanted variation (RUV), that adjusts for nuisance technical effects by performing factor analysis on suitable sets of control genes (e.g., ERCC spike-ins) or samples (e.g., replicate libraries).

RUVSeq essentially fits a generalized linear model (GLM) to the expression data, where your expression matrix $Y$ is a $m$ by $n$ matrix, where $m$ is the number of samples and $n$ the number of genes. The model boils down to

$Y = X*\beta + Z*\gamma + W*\alpha + \epsilon$

where $X$ describes the conditions of interest (e.g., treatment vs. control), $Z$ describes observed covariates (e.g., gender) and $W$ describes unobserved covariates (e.g., batch, temperature, lab). $\beta$, $\gamma$ and $\alpha$ are parameter matrices which record the contribution of $X$, $Z$ and $W$, and $\epsilon$ is random noise. For subset of carefully selected genes (e.g., ERCC spike-ins, housekeeping genes, or technical replicates) we can assume that $X$ and $Z$ are zero, and find $W$ - the "unwanted variation" in your sample.


fuse fastq files with multiple records I'd like to fuse ~50 pairs of corresponding fastq files on per sequence bases. I'd be happy for a solution as a Linux script, R script (if feasible with huge files), or any free software that runs preferably under Debian Linux (otherwise Windows). Aim is preprocessing for a Linux pipeline for NGS data. Ideally both the sequence and quality data are fused.

The aim of this request is the following: there is an existing pipeline for downstream mapping of immunoreceptor sequences (migec/mixcr/vdjtools etc.). Previously we got Illumina reads (read 1, read 2 = reverse), with the index and UMI sequences attached to one of the reads => can be alligned and processed to habe 1 biological useful sequence with UMI and index at one end. NOW we get read 1, read 2 and seperately "read 3" which is just the UMI+Index reads, and which was attached by the core facility to read 1 using software we do not have access to. Since the previous pipeline works well, we prefere to use it in the future, but are open for suggestions how to use the index reads with migec/mixcr/vdjtools otherwise.

Each of the 2 files contains the same number (millions) of short (~300 bp) and corresponding DNA reads:

File_A.fastq:

@record1 / aaa... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

File_B.fastq:

@record1 / ttt... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

The target file should look like:

@record1 / aaa...ttt... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

(of course the aaa...ttt... is just an example)

I was looking for a solution on Stackexchange and Google, but the solutions, I found do not work for me because

manual processing is not an option given the number of sequences
I need to fuse the files on a per sequence basis, i.e. a simple Linux command like cat file1 file2 is not what I need (as far as I understand, this would give me (record1_A record2_A.... record1B record2_B).
R, package shortread, readFastq + writeFastq would be perfect, but cannot handle such large numbers of input lines

Maybe I did not use the right search terms, but in any case I'd be happy for a solution

Based on the comments of the question I thought to put together this answer using GNU coreutils which should exist on pretty much all systems.

The idea is to paste the 2 files together with paste (assuming they are in the same order and each line correspond to the same line in the other file).

Since only the sequence and quality of both files should be kept, the lines starting with @ and + need to be made empty for one of the files. This can be achived with sed, e.g. sed 's/^@.*//' empties all lines starting with @ but keeps the empty line.

To modify the second file on the fly, one can using process substitution <().
Everything strung together should be something like this:

file1=File_A.fastq
file2=File_B.fastq
paste -d "" $file1 <(cat $file2 | sed 's/^@.*//' | sed 's/^+.*//') > target.fastq

I haven't tested it thoroughly, but could be a good starting point.

If you have a lot of files to paste, you might however want to process them first (not on the fly) using:

for file in $(ls *.fastq); do
  cat $file | sed 's/^@.*//' | sed 's/^+.*//' > ${file%.fastq}_cleaned.fastq
done

This creates File_A_mod.fastq from File_A.fastq etc.
Best to move them in a separate directory (make sure the first one is the original instead of the _cleaned one since it serves as the base).
You should have the following in the folder: File_A.fastq File_B_cleaned.fastq File_C_cleaned.fastq ....
Those can be than pasted with paste -d "" *.fastq > target.fastq.


Simple Affy Differential Expression I think I have a fairly easy-to-solve problem.

I have an expression table generated from an Affy set. It has one row of gene names (already symbols) and then 180 sample rows. 181 rows total.

So I made another excel file with one row with three different phenotypes I want to compare. That row has a label ("classes") and 180 names. 181 rows total.

In R, I can import each of these

dat <- read.csv("expression.csv")
dat <- data.matrix(dat)
pheno <- read.csv("pheno.csv")
pheno <- data.matrix(pheno)

I can then fit them

fit <- lmFit(dat,pheno)

And ask for an eBayes from limma

fit <- eBayes(fit)

Now, when I view the topTable,

table <- topTable(fit)

My first column has numbers (like "834","1142"), when I'd expect to have gene names there (like "ACTR", "CHRNA1").

What am I doing wrong? Or, rather, am I doing anything right?

Thanks!

What do the first few lines of your lmFit input variables look like (i.e. dat and pheno)? Are the rows numbered, or do they use the actual symbols?

By default, I don't think read.csv expects a file to have row numbers (even though write.csv puts them in by default). Try adding a row.names parameter to the expression matrix read call:

dat <- read.csv("expression.csv", row.names=1)


How to draw a 3d structure using 3d coordinates from pdb file? I am trying to create a 3d structure from 3d coordinated from a pdb file.
I am trying to draw 3d coordinates from pdb file "ATOM 1 N VAL A 14 22.692 -27.591 -39.402" and using webgl javascript.
I've managed to get the 3d coordinates from a pdb file but i don't know the exact logic that is used to create 3d structure from the coordinates found in pdb file.
I tried those coordinates and managed to create a cloud of 3d points but it is not anywhere near the structure that is found on pdb site.
So wondering how the structure is created on the pdb website?

JSmol sounds like what you want - it'll let you display 3D structures from PDB files onto webpages using JavaScript. If you just want to view the structure rather than display them on webpages, there's many programs that exist already like Swiss-PdbViewer and Pymol. There's even web-interfaces that let you upload a PDB file and they display it for you like Jena3D Viewer.


Question about the dots on Quartile groups in boxplot I have Microarray Normalized Expression data for a specific Gene. It looks like below in a dataframe B

SampleID    Gene    Type
Sample1     5.02    Tumor
Sample2     5.06    Tumor
Sample3     5.1     Tumor
Sample4     5.11    Tumor
Sample5    5.127    Normal
Sample6     5.12    Normal
Sample7    5.138    Normal
Sample8    5.149    Normal

I see that the minimum expression value in the table is 5.0 and maximum expression value is around 5.9.

I wanted to show the expression between two conditions with a boxplot and used following code.

q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               add = "jitter", ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

This gave a plot like this

But when I remove the jitter from the code,

q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

I see many black dots in the middle of the boxplot like this

May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?

You add the points with geom_point(). Just remove it and you will get your "empty" boxplot.

q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))

Unfortunately I couldn't use stat_compare_means(method = "t.test") and     stat_n_text(), but this two just add the labels.


Downloading an open-source MAF with tumor and normal BAMs? I'm writing an R package that manipulates MAF files. Normally, I generate a small example of a file and use this for testing. (e.g. a subsetted BAM or MAF).

In this case, I would need a MAF and a tumor/normal BAM pair. My problem is I don't know where I could download an open-source MAF with said BAMs.

Where could researchers download files like this, and not violate privacy laws (i.e. they must be publicly-available BAMs)?

There's the NCI Genomic Data Commons:

  For all submitted sequence data, including BAM alignment files, the GDC generates new alignments in BAM format using the latest human reference genome GRCh38 with standard alignment pipelines. Using these standard alignments, the GDC generates high level derived data, including normal and tumor variant and mutation calls in VCF and MAF formats, and gene and miRNA expression and splice junction quantification data in TSV formats.

But that needs an NIH account for individual-level data, which might be too restrictive for you.

There's also the Texas Cancer Research Biobank (7 human cancer patients with matched normal):

  These tumor and normal specimen sequence files (FASTQ and BAM), somatic variant calls (VCF and MAF) and germline MAF files are available for each consented patient as described on the TCRB Open Access Privacy page.

That also requires account registration, but the registration process is automated and much simpler (requires just a name, and acceptance of usage conditions). There's a published paper about that dataset.

I had some vague feeling that another repository existed (created within the last year), but I can't seem to find it from quick web searches.


Your alignment is

SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF

For each position in the alignment you calculate the score for that alignment. For position 1 we'd look up S vs R in the matrix and find a score of -1. You continue doing this until you hit the first -, which is not in the matrix. The first - is a gapopening, each consequent - in a series of -'s counts as a gap extension, instead of an opening.

---G--
IKAGDS

gives
-12 -2 -2 6 -12 -2


How to quickly and robustly convert between mmCIF and PDB? There is already a question on PDB/CIF to MMTF, however what is a robust way to programmatically go between PDB and CIF files?

For example I can use a python script from this gist that relies on Biopython.

id="4ckh"
wget http://www.ebi.ac.uk/pdbe/static/entry/download/$id-assembly-1.cif.gz
gunzip $id-assembly-1.cif.gz
python cif2pdb.py $id-assembly-1.cif $id.pdb

This process yields a pdb structure as shown below.

The issues are:

Requires python/Biopython to be installed.
The output isn't perfectly clean.
The gist is relatively undeveloped and hasn't got much evidence of testing.

Converting a PDB file to some version of an mmCIF file is always possible provided the PDB file doesn't have any "problems". However grouping identical molecules into the same entity is non-trivial so the Biopython and BioJulia mmCIF writers (which I wrote) do not try and do this, simply assigning new entites for each molecule. These files will also be missing other things like the metadata parts of the mmCIF dictionary not involving the atoms.

In Biopython:

from Bio.PDB import PDBParser, MMCIFIO
p = PDBParser()
struc = p.get_structure("", "file.pdb")
io = MMCIFIO()
io.set_structure(struc)
io.save("file.cif")

In BioJulia:

using BioStructures
struc = read("file.pdb", PDB)
writemmcif("file.cif", struc)

Going the other way is, however, not always possible. You can run out of atom serials or chain IDs to use in the output PDB file (this is one reason for the switch to mmCIF files in general). The gist you have linked is a way to get round this.

If this is not a problem for your structure you can invert the above snippets to get a decent PDB file - see the relevant docs. Note that in both cases it will use the auth_seq_id and auth_asym_id fields to get the residue number and chain ID rather than the relabelled label versions. The test suites of Biopython and BioJulia do test these conversions pretty thoroughly.

There is also pdb-tools if you do not want to install Biopython, though both involve a simple pip install so it's not much difference.

For some speed benchmarks of various software see here and here.


Compare alignment quality of multiple sequencing runs aligned against the same reference genome I have run Oxford Nanopore Technologies' MinION sequencing on the same DNA sample using three flowcells, each aligned against the same reference genome (E.coli K12 MG1655) using both BWA MEM and GraphMap and stored as BAM files.

How can I quantitatively and efficiently analyse the quality of alignment (percentage identity, insertion rate, deletion rate) of each of these files?

Qualimap will do this for you.

Go to qualimap.bioinfo.cipf.es
Run qualimap (default params are fine) on each BAM file
Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.

One thing to watch out for (you don't mention it in your question, but just in case) is that you cannot directly compare Q scores - these are a bit of a mess and calculated very differently in each piece of software.

Unsolicited suggestion: you might also try NGM-LR for mapping MinION data. We've found it beats the others for our data (though we map to a distant reference).


Can SBOL encode DNA Assembly linkers? Where in an SBOL file could you specify the linkers that would be placed as prefixes / suffixes between parts for assembly? Lets say you have a final construct of your design and want to assemble the component parts - could you encode what linkers you used between each component, ie as an annotation? Would this be different between SBOL versions?

Every version of SBOL allows you to specify prefixes/suffixes as annotations on a sequence.

In SBOL 3, the construct would be a Component of type DNA (https://identifiers.org/SBO:0000251), and the linker a SequenceFeature on that Component. The role of the sequence feature might depend on the particular assembly you are describing, and whether you're talking before or after assembly (in those cases where the process changes the sequence), but a good one in many cases would be Nucleotide Binding Site (https://identifiers.org/SO:0001655)
In SBOL 2, the representation is the same, except the construct would be a ComponentDefinition and the linker a SequenceAnnotation.
In SBOL 1, the representation is still the same, except the construct would be a DnaComponent and the linker an SBOL 1 SequenceAnnotation (which has different fields than in SBOL 2).

If the assembly method does not change the sequence, then all three can represent this pretty much equally well.  If you want to represent a sequence that is modified by assembly or the plasmid carriers from which parts are extracted, then SBOL 3 is the best because it cleanly represents deletions during assembly via the sourceLocation field.


How I know which gene is a good predictor in this neural network or not? I have 25 highly differentially expressed genes among and patients to chemotherapy.
I have made a neural network of these genes. The accuracy of the model is 0.73 but I don't know from 25 genes how I could select the best predictors or remove the weaker ones.
> # Convert column from /to 1/0
> = as.numeric(data$)-1
> data1 = cbind(,df)

> head(data1[,1:4])
   ALB CALML5 CALML6
A2       0 4.64   5.77   3.16
A3       0 7.12   9.13   5.13
A4       0 7.84   7.09   4.85
A6       0 7.90   9.64   6.52
A7       0 5.84   6.26   0.00
A8       0 4.83   5.77   2.15
>
> library(caTools)
> set.seed(101)
> # Create Split (any column is fine)
> split = sample.split(data1$, SplitRatio = 0.70)
> # Split based off of split Boolean Vector
> train = subset(data1, split == TRUE)
> test = subset(data1, split == FALSE)
> feats <- names(data[,-1])
> # Concatenate strings
> f <- paste(feats,collapse=' + ')
> f <- paste('TG ~',f)
> # Convert to formula
> f <- as.formula(f)
> f
~ ALB + CALML5 + CALML6 + CDK6 + CHGA + CREB3L3 + CSF3 +
    CXCL5 + CXCL6 + CXCL8 + DKK1 + EGLN2 + FGF20 + FGF4 + IL11 +
    IL1B + IL6 + ISG15 + KLK5 + KRT17 + MMP3 + MMP7 + PKP1 +
    PTGS2 + S100A7A
> #install.packages('neuralnet')
> library(neuralnet)
Warning message:
package ‘neuralnet’ was built under R version 3.5.2
> nn <- neuralnet(f,train,hidden=c(25,2,1),linear.output=FALSE)
> View(test)
> # Compute Predictions off Test Set
> predicted.nn.values <- compute(nn,test[2:26])
> # Check out net.result
> print(head(predicted.nn.values$net.result))
             [,1]
A4  0.26890819169
B5  0.94836576493
B7  0.19753554824
B8  0.99790838807
B9  0.99764782073
B12 0.05230351043
> predicted.nn.values$net.result <- sapply(predicted.nn.values$net.result,round,digits=0)
> table(test$Private,predicted.nn.values$net.result)

    0 1
  0 4 3
  1 1 7

> accuracy<-(table(test$Private,predicted.nn.values$net.result)[1,1]+table(test$Private,predicted.nn.values$net.result)[2,2])/(sum(table(test$Private,predicted.nn.values$net.result)))

> accuracy
[1] 0.7333333333
>

Yes, there is a formal solution. You want a heat map of the layers weighting and this is the information sought. It is more difficult to interpret than other supervised learning algorithms, but that is the nature of neural networks.

I know Sci-test Learn (Python), but I don't know R for neural networks, so I can't advise you on the code, only that it will be a few lines of code and a graphical (relative) output. I imagine Tensorflow will have a really brilliant heat map layer weighting.


Where can I find the genotype frequencies for the SNP's that are tested at 23andMe, Ancestry and FTDNA/MyHeritage? I want to add the frequencies of the genotypes to a baseline of raw DNA data I've got from the DNA testing companies 23andMe, Ancestry and MyHeritage (tested at all). All in build 37 human genome.
Where can I find that information? Checked HapMap project to no avail, can't align with the rsId or positions used in the text files.
Thanks in advance for your help

gnomad is probably the most comprehensive allele frequency database that exists at the moment.
Depending on how many variants you need to query, you can either manually input them into the search bar, or if you have a lot more, then you can also directly download the vcfs under the downloads tab.
The default for V2 gnomad set is on the gr37 build.


What are some good practices to follow during EPIC DNA methylation data analysis? I recently got some EPIC DNA methylation data and I was wondering what are some good practices to follow?
I am interested in knowing about normalization and differential analysis. Thank you.

EPIC data can be processed in the same manner as the previous iteration of methylation array data from Illumina (450k). This means that starting with .idat files, normalization should be performed (for example, via the minfi package). A recent paper from the creators of minfi is particularly helpful because it makes clear that normalized EPIC data from their package can be immediately compared against, for example, level 3 TCGA data.

After that, I suggest using the manifest to attach genomic coordinates to your probes and segregating them into functional regions. By testing differential methylation in only regulatory regions, for example, you can increase the statistical power by reducing the overall number of tests to the ones you expect to yield major differences.

There are existing packages out there for differential methylation analysis, but without knowing your replicate structure or aims, it is difficult to point you in the right direction.


Matching files by chromosome during wdl scatter gather I'm trying to implement a polygenic score pipeline in WDL and I'm quite new to this pipeline management.

The first step of the pipeline takes a file of GWAS summary statistics and splits it by chromosome.

The corresponding wdl task:

task split {
    input {
        File gwas
        String output_prefix
    }

    command {
        ./splitGwas -i ${gwas} -o ${output_prefix}
    }

    output {
        Array [File] gwas_by_chr = glob("${output_prefix}_*.assoc")
    }
}

The next part of the process is to compute posterior effects of SNPs by chromosome. The inputs to this task are going to be the files split by chromosome in the previous step and an LD matrix for the chromosome.

The corresponding wdl call would be something like this:

Array [File] ld_matrices = read_lines(file_of_ld_matrices_by_chr)

scatter(pair in zip(split.gwas_by_chr, ld_matrices)) {
    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = pair.left,
            ld_matrix = pair.right,
            output_prefix = out
    }
}

Since I'm grabbing the split files using glob in the previous step, I do not particularly know which element in the array corresponds to which chromosome. So I would be making a mistake if I zip the summary stats for one chromosome with the LD matrix corresponding to a different chromosome.

The programmatic way would be to iterate the scatter over chromosome name but then I would lose the ability to specify the dependency between the two tasks.

Is there a better way to do this?

Cross posted on biostars (https://www.biostars.org/p/434582/) so pardon me if that is discouraged.

You could use the sub function to obtain the chromosome from the globbed file names and then use that to retrieve the correct LD matrices from a Map (see read_map or read_json).

eg.:

Map[String, String] ld_matrices = read_json(file_of_ld_matrices_by_chr)

scatter(gwas_for_chr in split.gwas_by_chr) {
    String chr = sub(sub(gwas_for_chr, prefix, ""), "\.assoc$", "")

    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = gwas_for_chr,
            ld_matrix = ld_matrices[chr],
            output_prefix = out
    }
}

The file_of_ld_matrices_by_chr would have to look like this:

{
  "chr1": "/path/to/file1",
  "chr2": "/path/to/file2"
}


correlation between transposon and gene I have two matrix identical in N° of rows and columns. In the first matrix rows indicate transposable elements and columns  indicate samples. the matrix show the degree of expression of each element in samples:

M1 <-

            elements          ova       testes       retina   optic_lobe       suckers
1   LINES/3           1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   SINES             1.836210e+02   355.808323   577.807528   341.048786   99.902702
3  DNA TRANSP         9.181051e+00   178.661200  4387.148414    53.907711   353.413028

In Second matrix rows indicate genes and columns indicate samples. the matrix show the degree of expression of each gene in samples:

 M2 <-
            elements          ova       testes       retina   optic_lobe  suckers
1   OCBIM_22028242mg 1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   OCBIM_22008718mg 1.836210e+02   355.808323   577.807528   341.048786   99.902702
3   OCBIM_22009482mg 9.181051e+00   178.661200  4387.148414    53.907711   353.413028

I have to make an element-gene correlation. I have to check, for example, if the element A expressed much in the tissue X has the corresponding gene A expressed very, little, or nothing in the tissue X

fOR example I'm expecting this result:

M1,2 <-
                elements          ova       testes       retina   optic_lobe  suckers
    1   TE/ gene1                  0.3          1          0.5         -1         0.3
    2   TE/gene2                   0.1        0.6            0         -1         0.9
    3   TE/gene3                     1          1          0.8        -0.2       -0.3

Ok I think I roughly get what you want to get. Below I simulate something that's like your data because you provided two identical matrices.

#3 TE, 10 tissue, random values
# number of tissue
M1 = matrix(rnorm(30),3,10)
rownames(M1) = paste("TE",1:3,sep="_")
colnames(M1) = paste("Tissue",letters[1:10],sep="_")
## make gene matrix, , random values
M2 = matrix(rnorm(30),3,10)
rownames(M2) = paste("Gene",1:3,sep="_")
colnames(M2) = paste("Tissue",letters[1:10],sep="_")
## here we make TE1 highly correlated to Gene1
M2[1,] = M1[1,] + rnorm(10,0,0.3)
## here we make TE2 highly anti-correlated to Gene2
M2[2,] = -1*M1[2,] + rnorm(10,0,0.3)
#here we calculate for each TE, its correlation with each gene.
#credit to @haci
COR = apply(t(M1),t(M2))
# you can look at the correlation
COR

What you need is the last three rows. Remember, you can only calculate correlation between the gene and the TE, and you do this across all tissue. For example, if gene1 and TE1 has a high correlation, this tells you TE1 expression can predict gene1 expression in most tissue.

It's meaningless to calculate correlation between one data point and another. You need more than 1 value. For example, if you do cor(10,5), it doesn't make sense.

Hope this is what you are trying to do..


How to calculate average minor allele frequency (MAF) difference between two populations (i.e. European, Hispanic) The data is present in the PLINK format .bed (individual genotypes), .bim (genetic markers) and .fam (sample IDs and disease phenotype file). I have also separate text files containing information on subjects related to two populations. (i.e. European and Hispanic)

You need to report allele frequencies using plink for some European samples
Do the same for Hispanic samples
Use a software to subtract the two allele frequencies obtained above and obtain the result.


Error in cor(exprs(gse), use = "c") : no complete element pairs I'm trying to plot a heat map for a "GSE133399", I retrieved the GSE using GEOquery package and then assigned to an object and then tried to plot a heatmap, but I was not successful, I used the following code in RStudio:
install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)
## change my_id to be the dataset that you want.
my_id <- "GSE133399"
gse <- getGEO(my_id)
length(gse)
gse <- gse[[1]]
gse
pData(gse) ## print the sample information
fData(gse) ## print the gene annotation
exprs(gse) ## print the expression data

## exprs get the expression levels as a data frame and get the distribution
summary(exprs(gse))
exprs(gse) <- log2(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
library(dplyr)
sampleInfo <- pData(gse)
sampleInfo
## source_name_ch1 and characteristics_ch1.1 seem to contain factors we might need for the analysis. Let's pick just those columns

sampleInfo <- select(sampleInfo, source_name_ch1,characteristics_ch1.1)

## Optionally, rename to more convenient column names
sampleInfo <- rename(sampleInfo,group = source_name_ch1, patient=characteristics_ch1.1)
sampleInfo
install.packages("pheatmap")
library(pheatmap)
## argument use="c" stops an error if there are any missing data points

corMatrix <- cor(exprs(gse),use="c")

I get this error:
Error in cor(exprs(gse), use = "c") : no complete element pairs

I also  tried this:
corMatrix <- cor(exprs(gse),use="c", use="complete.obs")

I get this error:
Error in cor(exprs(gse), use = "c", use = "complete.obs") :
  formal argument "use" matched by multiple actual arguments

I want to achieve this plot:

The GSE133399 dataset contains both ATAC-seq and RNA-seq data, none of which can be queried with GEOquery which works only for microarray data.
exprs(gse) ## print the expression data returns an empty object as you might have noticed, therefore the errors in the lines that follow this command.
On the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133399 at the Supplementary file section you can download FPKM-transformed count data and use these. I would be careful though as FPKM is not ideal for inter-sample comparison. I personally always (if feasable) download the raw data and then process myself, but for a guess on what is going on the FPKMs probably will do.


What is the `table` parameter in the UCSU ldHgGene tool? I want to lift over the annotation from E. coli release 29 to E. coli release 42. I have followed the steps for an annotation liftover mentioned here. I got the chain file for the final step and ran:

liftOver -gff Escherichia_coli.GCA_000967155.1.gtf Chromosome.chain conversions.gtf unmapped.gtf

However, I get the following warning:

WARNING: -gff is not recommended.
Use 'ldHgGene -out=<file.gp>' and then 'liftOver -genePred <file.gp>'

I would like to heed the warning. However, I don't know how to use the ldHgGene tool. Here is its help page:

ldHgGene - load database with gene predictions from a gff file.
usage:
     ldHgGene database table file(s).gff
options:
     -bin         Add bin column (now the default)
     -nobin       don't add binning (you probably don't want this)
     -exon=type   Sets type field for exons to specific value
     -oldTable    Don't overwrite what's already in table
     -noncoding   Forces whole prediction to be UTR
     -gtf         input is GTF, stop codon is not in CDS
     -predTab     input is already in genePredTab format
     -requireCDS  discard genes that don't have CDS annotation
     -out=gpfile  write output, in genePred format, instead of loading
                  table. Database is ignored.
     -genePredExt create a extended genePred, including frame
                  information and gene name
     -impliedStopAfterCds - implied stop codon in GFF/GTF after CDS

From what I understand, I am trying to give the liftOver tool the annotation in GenePred format instead of GTF format. So I need to convert the GTF format into a GenePred format and that's what the ldHgGene tool is for. So the command should be something like:

ldHgGene Escherichia_coli.GCA_000967155.1.gtf -out=Escherichia_coli.GCA_000967155.1.gp

but this returns nothing. The -out option does mean that the database parameter is not needed but the tool still seems to need the parameter table.

What is the table parameter?

I stumbled across this issue after getting the same warning message from liftOver. By using the -out option of lgHgGene, the database parameter is ignored but must still be provided. The table parameter should be the second column from your GTF file. Note that you must also specify the -gtf option (otherwise it expects GFF format). In my case, the GTF file looks like:

chr1    mm10_knownGene  exon    3073253 3074322 0.000000        +       .       gene_id "ENSMUST00000193812.1"; transcript_id "ENSMUST00000193812.1";
chr1    mm10_knownGene  exon    3102016 3102125 0.000000        +       .       gene_id "ENSMUST00000082908.1"; transcript_id "ENSMUST00000082908.1";
chr1    mm10_knownGene  exon    3206523 3207317 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3213439 3215632 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3205901 3207317 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
chr1    mm10_knownGene  exon    3213609 3216344 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
...

By default, ldHgGene outputs a leading 'bin' column which won't work with liftOver; this can be avoided with the -nobin option. So my command was:

ldHgGene -gtf -nobin -out=mm10.GENCODE_VM20.gp ignored mm10_knownGene mm10.GENCODE_VM20.gtf

The resulting genePred output can be fed directly into liftOver -genePred:

ENSMUST00000193812.1    chr1    +       3073252 3074322 3074322 3074322 1       3073252,        3074322,
ENSMUST00000082908.1    chr1    +       3102015 3102125 3102125 3102125 1       3102015,        3102125,
ENSMUST00000162897.1    chr1    -       3205900 3216344 3216344 3216344 2       3205900,3213608,        3207317,3216344,
ENSMUST00000159265.1    chr1    -       3206522 3215632 3215632 3215632 2       3206522,3213438,        3207317,3215632,


How BWA generate index files? I need to know how BWA generate bw and sa in less memory usage ? Do they save all rotations temporary ? I need to know do BWA stores all rotations so that it can sort the rotations alphabetically ?

Thanks.

Papers

The space-efficient BWT construction algorithm is in bwt_gen.c. One of the first things you see there is the copyright line:

  Copyright (C) 2004, Wong Chi Kwong.

With that name, you can find the relevant papers:

Lam et al: Compressed indexing and local alignment of DNA (Bioinformatics, 2008). This paper describes the BWT-SW alignment algorithm.
Hon et al: Practical Aspects of Compressed Suffix Arrays and FM-index in Searching DNA Sequences (ALENEX 2004). This was the first paper to demonstrate that it's feasible to build the FM-index for a human genome.
Hon et al: A Space and Time Efficient Algorithm for Constructing Compressed Suffix Arrays (Algorithmica, 2007). This paper contains the theoretical description of the construction algorithm.

The third paper only mentions BWT briefly in the last section. It deals with compressed suffix arrays, which are space-efficient representations of the inverse function of LF-mapping. But if you have a representation of (the inverse of) the LF-mapping, you have a representation of the BWT. Hence the algorithm can be adapted for BWT construction.

Algorithm

When building BWT, you usually append a unique character $ to the end of the text, either implicitly or explicitly. This way you can reduce sorting rotations to sorting suffixes. You sort the suffixes in lexicographic order, and then you output the preceding character for each suffix in that order to form the BWT.

Assume that you already have the BWT for string X, and you want to transform it into the BWT of string cX, where c is a single character. You can do that in three steps (this assumes 0-based array indexing):

Find the position i of the full suffix X, and replace the $ in the BWT with the inserted character c (set BWT[i] = c).
Determine the number of suffixes j that are lexicographically smaller than cX with j = LF(i, c) = C[c] + BWT.rank(i, c). Here C[c] is the number of suffixes that start with a character smaller than c, and BWT.rank(i, c) is the number of occurrences of c in the prefix BWT[0..i-1].
Insert $ between BWT[j-1] and BWT[j].

By iterating this, you can build the BWT for any string slowly but space-efficiently. The algorithm of Hon et al. is based on the same idea, but they make it faster by inserting multiple characters at the same time. Given the BWT of string X and another string Y, they produce the BWT of the concatenation YX. The details are quite complicated, but you can find them in the third paper and in the source code.


How I can change the name of multiple files at once in R or terminal? I have 200 .vcf files in a folder with long names like

LP6005409-DNA_E03_vs_LP6005408-DNA_E03.snp.pass.vcf

How can I change the name of each file to, for example A1.vcf ....A200.vcf?

Renaming these files manually takes ages, help please.

> file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))
Error in file.rename(list.files(pattern = "vcf"), paste0("A", 1:length(list.files(pattern = ".vcf")),  :
  'from' and 'to' are of different lengths

[fi1d18@cyan01 new]$ for i in *.vcf; do echo mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
mv LP2000325-DNA_A01_vs_LP2000317-DNA_A01.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
mv LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
[fi1d18@cyan01 new]$ for i in *.vcf; do mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
(standard_in) 2: syntax error
mv: try to overwrite `A.vcf', overriding mode 0500 (r-x------)?

Use perl-rename. This is usually called rename on Debian-based systems like Ubuntu or Mint, and perl-rename or prename on others. Assuming you have it as rename, you can simply do:

rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf

The -n causes rename to only print what it would do, without doing anything. So once you make sure it does what you want, run the command again without the -n.


Problem installing RATT I have followed the steps here to install RATT, i.e.:

svn co "https://svn.code.sf.net/p/ratt/code/" ratt-code

However, when I check my installation by running '$RATT_HOME/start.ratt.sh', it says:

zsh: no such file or directory: /start.ratt.sh

which makes sense once I check that the RATT_HOME variable is empty.

I downloaded their tuberculosis example, went into that directory and ran:

~/path/to/my/ratt/ratt-code/start.ratt.sh embl F11.fasta F11 SimStrain

and got the following output:

Please set the RATT_HOME variable.
At Sanger for bash it is RATT_HOME=/nfs/users/nfs_t/tdo/Bin/ratt; export RATT_HOME
At Sanger for tcsh setenv RATT_HOME /nfs/users/nfs_t/tdo/Bin/

What should I set my RATT_HOME variable to? The message only says what it should be set to at the Sanger Institute.

Usually the home of a program is where the executable is. In order to make the system aware you need to include it as a variable when the terminal is run.

So you need to add it to the file at ~/.bashrc of your computer (if you are using bash as shell [you can check it with echo $SHELL]). You need to add :

export RATT_HOME=path/where/RATT/is/installed/or/downloaded

As terdon says in a comment below:

  Note that echo $SHELL just prints the default login shell of the
  user, and not necessarily the shell that is currently running.  Also,
  ~/.bashrc is not read by macOS bash which starts login shells by
  default. In any case, this sort of global variable should be in
  ~/.profile, not .bashrc. For more of this sort of detail see
  Sequence of scripts sourced upon login and How to correctly add
  a path to PATH.


Why do some assemblers require an odd-length kmer for the construction of de Bruijn graphs? Why do some assemblers like SOAPdenovo2 or Velvet require an odd-length k-mer size for the construction of de Bruijn graph, while some other assemblers like ABySS are fine with even-length k-mers?

From the manual of Velvet:

  it must be an odd number, to avoid palindromes. If you put in an even
  number, Velvet will just decrement it and proceed.

the palindromes in biology are defined as reverse complementary sequences. The problem of palindromes is explained in this review:

  Palindromes induce paths that fold back on themselves. At least one
  assembler avoids these elegantly; Velvet requires K, the length of a
  K-mer, to be odd. An odd-size K-mer cannot match its reverse
  complement.

It is possible to construct graph with palindromes, but then the interpretation will be harder. Allowing only graphs of odd k-mers is just an elegant way to avoid writing a code for interpretation of a more complicated graph.


Can someone recommend software for reference-based genotype-estimation? Problem: I have low-coverage sequencing data from humans in vcf format. For each sample, I have PL scores / preliminary genotype calls from the software ATLAS (it's a bit like GATK for those who don't know). Because the samples are low coverage (mean 1x, with lots of missing data), I am looking to enhance the genotype calls by using a reference panel and genetic map, before go on to phase the samples. Briefly, the reference panel and genetic map can be used as a kind of prior to inform one about genotypes when there is little data present.

What I've done so far: So far I have been using Beagle4 with a command:

beagle4 gl=lowCoverageSamples.vcf ref=highCoverageReference.vcf map=geneticMap out=lowCoverageSamples.newCalls.vcf

However, Beagle4 is extremely slow, very memory inefficient and quite old. There have been many advances in phasing methods, which are now much more accurate, faster and memory efficient than the beagle4 model. However, I haven't come across any similar improved methods for the genotype re-estimation step.

Does anyone knew of any newer pieces of software which do such a thing, which may more faster and more accurate than Beagle4. It needs to have 2 things in particular:

1) Accept data input in PL score format, as to allow the appropriate representation of uncertainty in low coverage samples.

2) Produce a genotype call / genotype probability for each SNP in each individual as output.

In case anyone was interested, a new piece of software for rapidly imputing and phasing low coverage sequence data from genotype likelihoods has just been released. It is called GLIMSPE and is exactly what I was after.

It looks to be much faster (1200x speedup) compared to Beagle 4.0, and is more accurate at imputing low frequency variants. It is also more amenable to use with larger reference panels (such as the Haplotype Reference Consortium panel).


Split FASTQ and matching BAM into matching chunks I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).

How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?

My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.

EDIT: Here's some pseudocode of exactly what I'm trying to do:

# input: in.bam, in.fastq, chunk_size
i <- 0
for fastq_read in in.fastq:
    bam_read <- extract fastq_read.read_name from in.bam
    n <- i modulo chunk_size
    write fastq_read to out.n.fastq
    write bam_read to out.n.bam
    i <- i + 1

I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.

Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.

Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:

import pysam
import itertools

def get_read_id_fastq(ref_path):
    """Extracts the read ids from a fastq file."""
    read_ids = set()
    with open(ref_path, 'r') as ref:
        for line in ref:
            if line.startswith('@'):  # i.e if line is header
                # split the line on spaces, take the first element, remove @
                read_id = line.split(' ')[0].replace('@', '')
                read_ids.add(read_id)

    return read_ids

def get_read_id_bam(ref_path):
    """Extract the read ids from a BAM file."""

    read_ids = set()
    with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
        for read in ref:
            # query_name is the query template name
            read_ids.add(read.query_name)

    return read_ids

fastq_ids = get_read_id_fastq(fastq_path)
bam_ids = get_read_id_bam(bam_path)

Then take the intersection of these two sets and that's your common read ids.

common_ids = fastq_ids & bam_ids

The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively

chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
write_idx = {}

The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.

for line in file:
    read_id = # get line read_id
    if read_id in common_ids:
        chunk_num = chunk_cycle.next()
        write_idx[read_id] = chunk_num
        file_to_write_to = 'out.{}.bam'.format(chunk_num)
        # open this file or write to it if already open

When you go to do the next file when you find a read id is in the common set, you look up it's value in write_idx and this will give you the chunk number to write to.

The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).

Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.


bedtool intersect to compare .vcf file to reference genome I have a .vcf file generated using samtools. I want to know mutations in specific genes by comparing .vcf files with a reference genome. I cannot use the Ensembl VEP tool as my reference genome is not listed there. Can you please help?

If you want to filter a VCF file, you might want to give VCFtools a try.

Here is the documentation : VCFtools Manual

I am not a specialist of this tool but you could do as the following. First build a bedfile with the chromosome names and positions of your genes (start and stop) such as this :

chr   start   stop   (name)   (score)   (strand)

Then it should work with this command line :

vcftools --vcf file.vcf --bed positions.bed --out filtered.vcf

The output will be a VCF filtered on the provided regions


Correlating gene expression with qualitative variables I have a gene expression dataset that I want to investigate. Particularly, I would like to understand whether there is any correlation between each gene's expression and some quantitative or qualtitative data (say, correlation between gene 'XPTO' , body mass index, and race).

One possible way to test this would be through logistic regression, but is this a good approach or are there caveats that I should know about using such a statistic?

My question is the following: which methods would you advise to measure such correlations, and why?

(this post was crossposted on Biostars)

Logistic regression would generally be a bad choice for non-binary outcomes. In such cases, linear regression (or a GLM more generally) still works fine. You can already do that in the standard R packages for RNAseq (DESeq2, edgeR, and limma), where the fold-change is in whatever units you're measuring your quantitative trait in.


Allele count is zero I ran freebayes on bacterial data and got a vcf file. Some of the SNPs found have AC=0, but in the bam file I can see that there are alternative alleles, so in my understanding it supposed to be AC>0.

example:
in vcf file:

CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  55509   55508   55506   55507
scaffold1       78      .       A       C       9.27683e-13     .       AB=0;ABP=0;AC=0;AF=0;AN=3;AO=466;CIGAR=1X;DP=2254;DPB=2254;DPRA=0;EPP=20.9226;EPPR=29.3274;GTI=0;LEN=1;MEANALT=2;MQM=50.3412;MQMR=45.5379;NS=3;NUMALT=1;ODDS=2080.75;PAIRED=0.899142;PAIREDR=0.533932;PAO=0;PQA=0;PQR=0;PRO=0;QA=16347;QR=63107;RO=1783;RPL=203;RPP=19.7856;RPPR=291.853;RPR=263;RUN=1;SAF=465;SAP=1006.25;SAR=1;SRF=894;SRP=3.04075;SRR=889;TYPE=snp       GT:GQ:DP:AD:RO:QR:AO:QA:GL      .       0:126.704:776:626,149:626:22118:149:5255:0,-1217.7      0:126.704:834:663,170:663:23671:170:6004:0,-1224.79     0:126.704:644:494,147:494:17318:147:5088:0,-770.853

but with samtools tview:

Can someone tell me if this makes sense and why it happens?

The AC is the number of alternate alleles from called genotypes, you can see this in the vcf header:

##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of
alternate alleles in called genotypes

Assuming diploid, if we have 3 individuals, 0/0 , 0/1, 0/0, AC will be 1, because it detects 1 count of alternate allele.

If you look at your genotype calls, all are 0, so AC = 0


Have DNA motifs 6-12bp long, trying to get conservation scores I have about 200 short nucleotide motifs (6-12 bp in length) from the human genome, and I'm trying to see how conserved they are across vertebrates.

I was thinking that I'd need to make a bed file for each motif that lists all of its occurrences in the human genome. From there, I could map the beds to a bigwig files of PhastCons scores (essentially doing the reverse of what the PhastCons software was designed to do). Does that sound like the best approach?

I'm getting stuck at the step of going from motifs to bed files. I've tried using BLAST to find all occurrences of motifs, but their short length is causing issues.
I've tried messing with the e-value threshold, word size, and filter parameters, but I still don't get any hits.

Is there a work-around for this issue, or should I just rethink my entire approach?

In case you have only ACGT in your motifs

The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

In case you have IUPAC bases (N, W, etc.) in your motifs

I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.


Creating a count table using R I’m trying to develop a table from a series of lists generated using the Immunarch library to process TCR sequencing results. Each list is comprised of CDR3.aa (clone) information which are character strings and their count in a particular sample. The clones are short but vary between 7 and 20 (or more).  Each list has a header that identifies the sample sequenced. I have a list of 66 samples. Each sample can have several thousand clone strings. Not every clone is contained in every sample so the number of clones listed in the samples varies.  Here’s an example of the structure of a single sample list.

$PreAg_18_2
CDR3.aa           Count
CASSYGTAYTGELFF   1623
CASSRGDSDNSPLHF   1440
CASSREKAFF        1161
CSGMGALAKNIQYF     949
CSAYTGLSYEQYF      813
CASSLSLAVNSPLHF    634
CAIRDTPGSPQHF      574
CATGQVNTEAFF       555
CASSLKGQGGSPLHF    499
CASSYSRSPQPQHF     478

I want to combine the results in a single table showing the clone counts with all the clone strings listed on the y-axis and the sample ID listed on the x-axis. For example:

                10_pep_10_1     preAg_10_2      Dec_2_18_1  …...
CASSYGTAYTGELFF    1623         234             0
CASSRGDSDNSPLHF    1440         522             28
CASSREKAFF         1161         445             50
CSGMGALAKNIQYF      949         24              0
CASSYSRSPQPQHF      478         0               398
.
.

Currently, I'm working with R because that's what Immunarch is in, but I wonder if a) There's a better way to do this in Python, or b) If it's doable in r, then how do I generate the matrix matches.  The main problem is that not all the samples have the same number of clones, so I first have to separate each sample data before I can convert it to a data frame, and then how can I match the actual count with the actual clone.  I can't just merge the data frames.  Any suggestion is greatly appreciated.

Here's what i have for code right now.
Code:
    library(immunarch)
    library(stringr)
    library(plyr)

immdata = repLoad("/mnt/data/Development/Analysis_Scripts/input_files/")

all <- immdata$data

# Get list headers (names) and convert to df
sample.id <- names(all)
sample.id <- data.frame(sample.id)

# Get list of clones and filter for unique clones per list.
for (i in 1:length(all)){
    all[[i]]$Sample.ID<-names(all)[i]
    all[[i]]<-all[[i]][,c("CDR3.aa", "Clones")]
}

# all is the variable that contains the list of samples and clones.
all.split <- split(all, )

# make vector of all clones
all.clones <- unlist(all, use.names=FALSE)

# Removes clone repeats
all.clones.u <- unique(all.clones)

# convert list of clones to data frame
all.clones.u <- data.frame(all.clones.u)

At this point I have both the list of sample ids and the list of clones.  My problem is now matching each count with respective sample and  clone in the table.  Any suggestions?

You could give a try with dplyr's full_join() and purrr's reduce(). full_join() will make sure that you won't lose any data and reduce() will merge tables in your list from left to right, left being the first table of the list in the first iteration but would become merged table1+table2 in the second iteration and so on. This left to right "reduction" will also be helpful when renaming the columns of the resulting table.

Without your workspace I have attempted to replicate your data: There is a list of 3 data frames, each with 2 columns:

library(dplyr)
library(purrr)

df1 <- data.frame(CDR3.aa = c("A", "B", "C"),
                  Count = c(10,20,100))

df2 <- data.frame(CDR3.aa = c("A", "B", "D"),
                  Count = c(100,200,300))

df3 <- data.frame(CDR3.aa = c("D", "E", "F"),
                  Count = c(10,20,100))

my_list <- list(df1 = df1, df2 = df2, df3 = df3)

> my_list
$df1
  CDR3.aa Count
1       A    10
2       B    20
3       C   100

$df2
  CDR3.aa Count
1       A   100
2       B   200
3       D   300

$df3
  CDR3.aa Count
1       D    10
2       E    20
3       F   100

# "by" argument specifies which column to use when merging
my_table <- reduce(my_list, full_join, by = "CDR3.aa")

names(my_table) <- c("CDR3.aa", names(my_list))

> my_table
  CDR3.aa df1 df2 df3
1       A  10 100  NA
2       B  20 200  NA
3       C 100  NA  NA
4       D  NA 300  10
5       E  NA  NA  20
6       F  NA  NA 100

Clonotypes that are not present in table will be "NA".


Download multiple fastq files using fastq-dump I want to download the following fastq files at the same time in Salmon:
 - SRR10611214
 - SRR10611215
 - SRR10611215
 - SRR10611216
 - SRR10611217

Is there a way to do this using a bash for loop or fastq-dump? or prefetch

A sample code is given in the salmon documentation as follows. Source
#!/bin/bash
mkdir data
cd data
for i in `seq 25 40`;
do
  mkdir DRR0161${i};
  cd DRR0161${i};
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz;
  cd ..;
done
cd ..

This could be modified as follows.
#!/bin/bash
mkdir data
cd data
for i in `seq 14 17`;
do
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_2.fastq.gz;
done
cd ..

You can save the code as a shell script and run it from the linux terminal. For example bash download_fastq.sh


Mapping drug names to ATC codes I'm interested working with the medication information provided by the UK Biobank.  In order to get these into a usable form I would like to map them to ATC codes.  Since many of the drugs listed in the data showcase include dosage information, doing an exact string match between drug names is not very effective.  I've considered using something like fuzzywuzzy to do string matching between the medications in the data showcase and the ATC drug names but validating the matches could still be a laborious process.  Does anyone know of a tool that can match drug names to ATC codes or some other drug ontology?  If not, maybe there's a better way to do it that I haven't thought of.

The CART tool let's you upload a set of names and map them (optionally in a fuzzy way) to STITCH 4 identifiers, and then use those to map to ATC codes (using the chemicals sources download file). It's a bit indirect, and I'm not sure what CART will do with the dosage info you mention.


What is the ICGC normalized_read_count? I downloaded gene expression data (exp_seq) from the ICGC file browser.

For each sample and gene, the file contains a normalized_read_count.

What is that value? I couldn't find any information on the ICGC website. The values are definitly too low for TPM.

By reading this thread on seqanswers and by comparing the data to TCGA, I figured out

raw_read_count is the read count which you use as input for e.g. DESeq2. It has been estimated using RSEM
normalized_read_count is equivalent to the scaled_estimate from TCGA. This is the estimated fraction of transcripts made up by a given gene, as estimated by RSEM. Multiplying this value with 1e6 yields the TPM.


Accidental mapping of eukaryotic reads in a metagenomic dataset This is a question from /u/wipeyourmit on reddit. The original post can be found here.

  If I have a metagenomic dataset that contains reads from both
  eukaryotes and prokaryotes and then I annotate by running DIAMOND or
  HMMER against a bacterial database how much of a risk do I run of
  eukaryotic reads being annotated in the process?

  I was hoping to use the eggNOG mapper to search against the bacterial
  and archaeal databases and to exclude the eukaryotic portion of my
  dataset. Is the eukaryotic filtering something that I would need to do
  in a step prior to this?

If you know the eukaryotic contaminant present you could use bbsplit.sh from BBMap suite to split/bin reads first using a reference for that contaminant (into one or as many bins as reference sequences you provide).


Which is the amino acids classification when analyzing the solvent accessibility of a protein? I came across an article which has defined 3 categories: Buried [A, L, F, C, G, I, V, W], Exposed [ R, K, Q, E, N, D] and Intermediate [M, S, P, T, H, Y] when analyzing the solvent accessibility of a protein. But I could not find any reference validating the claim.

Another article has come up with 4 divisions:

B=completely buried (0-4% exposed)
b=partly buried (4-25% exposed)
e=partly exposed (25-50% exposed)
E=completely exposed (50+% exposed)

Which is authentic when dealing with solvent accessibility, and which should be used? Also for the second article [B,b,e,E], which amino acid corresponds to which class?

There is some confusion in the question.

Solvent accessibility of a residue, i.e. the area accessible to the solvent (water), depends on the 3D structure of the protein. It's not a property of  an amino-acid (chemical compound such as alanine, leucine, ...).

There may be a correlation -- amino-acids with hydrophobic side chains are more likely to be buried, but in general amino-acids (chemical compounds) cannot be classified as buried or exposed.

I don't have access to the first article, but I tried running the prediction server described in the second linked article, and I got a result

GSGMQIFVKTFRDRLRNLLPQTITLEVEPSDTIENVKAKIQDKEGIPPDQQVLIFSRKRL
CCCCEEEEEECEEEEEECCCCEEEEEECCCCCHHHHHHHHHHHHCCCHHHEEEEECCEEC
EEEBeBbBebbEbEbEebEEeEbeBEBeeEeeBeeBbEeBeEeEEbEeEeBeBbbEEeeb

...

Porter (Secondary Structure):
H = Helix   (DSSP classes H, G and I)
E = Strand  (DSSP classes E and B)
C = Coil    (DSSP classes S, T and .)

PaleAle (Relative Solvent Accessibility):
B = very buried      (<=4% accessible)
b = somewhat buried  (>4% and <=25% accessible)
e = somewhat exposed (>25% and <=50% accessible)
E = very exposed     (>50% accessible)

In one place in the sequence you can see QQ, predicted as eB.
So clearly glutamine is NOT assigned to a particular RSA class.

Dividing RSA (relative solvent accessibility, 0-100%) into categories is arbitrary, I don't think there is any widely used convention. The easiest way would be to just use numeric values.


How to append numbers only on duplicates sequence names? I have a reference database with contains 100s of sequences in fasta format. Some of these sequences have duplicate names like so:

>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName
atgc
>1_duplicateName
atgc

Is is possible to run through a large file like this and change the names of only the duplicates?

>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName_1
atgc
>1_duplicateName_2
atgc

And here’s a solution using R (with the Bioconductor):

fa = ShortRead::readFasta(infile)
ids = as.character(ShortRead::id(fa))
fa@id = Biostrings::BStringSet(make.unique(ids, sep = '_'))
ShortRead::writeFasta(fa, outfile)


Use conditional in workflow in Nextflow dsl2 How should I do:
...
if ($(params.aligner) == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
}
else {
    align_bwa(get_fq_info.out.fq_info)
}
...

Is this inside the main.nf file or inside the process.nf file? Anyway, the param.anything variable should already be available in the nextflow context when executed, thus, you could simply do:
if (params.aligner == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
} else {
    align_bwa(get_fq_info.out.fq_info)
}

Check out this main.nf of a nextflow pipeline of mine, it contains some working if/else snippets.

Here it is another StackOverflow question on conditional nextflow pipelines that might help.

Hope it helps.