text:
How do population genetics people define a population?
How do population genetics people define a population?
Do they define it as a layman will do? say Africans, Americans, or so? Or is there a more scientific way of doing so? For example, I think defining one's population as one's allele frequencies makes more sense, but has this point been discussed previously?

Do you talk about a general definition or you want to what what defines human populations?

It really depends on the study and the ethnicity. Mostly they use self-defined ethnicity and often they will require that three or four grandparents are of that ethnicity. Here's some detailed info on how the ethnicity was defined for the 1000 Genomes populations. You'll see that some populations (such as ASW) require all four grandparents to identify as that ethnicity, but some (such as ACB) only need three, and the cultural sensitivities in Japan meant that JPT participants weren't questioned directly but told that the study was only interested in people with Japanese grandparents. The mixed American super-population (AMR) is an exception to this rule due to the extreme admixture of ethnicity in these regions, and is based more in geography than ethnicity – this means that the individuals had to be from the place and have grandparents from the place.

A population is defined as a group of people (or organisms of the same species) living together in a geographical region. However, for human genetic studies, it is wise to use ethnicity as well (so group of people with the same ethnic background). See here how they define population in the 1000 genomes project.


text:
What is Allelic Imbalance
Can anyone help me explain what allelic imbalance is, hopefully shortly?
Surprisingly, we cannot find any introduction online.

Some more context would help. Where did you see the term? Is it with respect to a phenotype? A gene? A population? An individual?

You may also consider looking into allele-specific expression. That may turn up more information on allelic imbalance.

It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions).
There are a few blog posts here with some more information, and plenty of papers (google "Allelic imbalance") if you want to dig deeper.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Finding annotated counterpart after BLASTn with efetch (Biopython)
I am creating a pipeline for the identification of unknown transcripts. After a local BLASTn search of the transcripts, I have a large list of the respective hits with different genomes. I have the accession code of the respective genome, start/end coordinates, and the strand where the hit is found. Now, I would like to find if the genome where the hit was found is annotated within the coordinates returned by BLAST.
I have found the possibility of using efetch to search for each genome via Biopython, but this only returns the whole annotation, even though I limit the search to the specific coordinates with:
Entrez.efetch(db="nuccore", id="CP054847.1", strand="2", seq_start="326053", seq_end="326786", retmode="xml")

Of course, I could parse the XML and look for possible annotations within the given range, but not only would this be quite computationally intense since I expect around 10.000+ hits to be searched, but it also seems to return the same results for the plus/minus strand, even though the annotations should be strand-specific.
I have also intended to use the standalone tool of efetch, but it would be easier to integrate the search to the pipeline using Biopython.
Do you have any ideas why Entrez.efetch is returning the whole xml or is there any fast way of filtering the output? Do you have maybe another alternative/idea for this step? I would appreciate any suggestion!

can you find any difference between https://www.ncbi.nlm.nih.gov/nuccore/CP054847.1/ and https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP054847.1 database entries ?

Lets assume you can get the blastn output in a BED format:
genomeA_contig001    hit_start    hit_end    query_transcript    score   strand
genomeB_contig021    hit_start    hit_end    query_transcript    score   strand

which will contain all the unknown hits to all genomes.
All you need is a set of GTF files specific for genomeA, genomeB and so on where you add the suitable prefix (genomeA_ for GTF describing genes/exons in that genome and so on). For transcript hits you may start with just the annotated exons. Next run:
bedtools intersect -wao -a all.transctipts.hits.bed -b all_annotated_exons.gtf my_result.tsv
to get all the overlaps with number of overlapping bases.
Filter and sort the my_result.tsv possibly importing it into a data frame using pandas or polars.

It seems like the seq_start and seq_end are sensitive to the direction. Try swapping the two.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Snakemake and SLURM -unable to use it with with --cluster
I'm trying to automate some of my previous bash scripts with snakemake. Below is one of the recipe I created:
import os
shell.executable("/bin/bash")

workdir: "/proj/uppstore2018136/publication/jupyterNotebooks"

ENSEMBL_RELEASE='94'
ENSEMBL_GRCm38_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/fasta/mus_musculus/dna'
ENSEMBL_GRCm38_GTF_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/gtf/mus_musculus'
GTF_FILE='Mus_musculus.GRCm38.'+ENSEMBL_RELEASE+'.gtf'

rule all:
    input:
        "genome/genome.fa",
        "genome/"+GTF_FILE,
        "indexes/bowtie2/reference.1.bt2",
        "indexes/bowtie2/reference.2.bt2",
        "indexes/bowtie2/reference.3.bt2",
        "indexes/bowtie2/reference.4.bt2",
        "indexes/bowtie2/reference.rev1.bt2",
        "indexes/bowtie2/reference.rev2.bt2",
        "indexes/hisat2/reference.1.ht2",
        "indexes/hisat2/reference.2.ht2",
        "indexes/hisat2/reference.3.ht2",
        "indexes/hisat2/reference.4.ht2",
        "indexes/hisat2/reference.5.ht2",
        "indexes/hisat2/reference.6.ht2",
        "indexes/hisat2/reference.7.ht2",
        "indexes/hisat2/reference.8.ht2"

rule getGenome:
    output:
        output1="genome/genome.fa",
        output2="genome/"+GTF_FILE
    shell:
           '''
           curl -o genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz {ENSEMBL_GRCm38_BASE}/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           gunzip genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           mv genome/Mus_musculus.GRCm38.dna.primary_assembly.fa genome/genome.fa
           curl -o genome/{GTF_FILE}.gz {ENSEMBL_GRCm38_GTF_BASE}/{GTF_FILE}.gz
           gunzip genome/{GTF_FILE}.gz
           '''

rule bowtie2Build:
    input: "genome/genome.fa"
    params:
        basename="indexes/bowtie2/reference"
    output:
        output1="indexes/bowtie2/reference.1.bt2",
        output2="indexes/bowtie2/reference.2.bt2",
        output3="indexes/bowtie2/reference.3.bt2",
        output4="indexes/bowtie2/reference.4.bt2",
        outputrev1="indexes/bowtie2/reference.rev1.bt2",
        outputrev2="indexes/bowtie2/reference.rev2.bt2"
    threads: 20
    shell:
        '''
        module load bowtie2
        bowtie2-build --threads {threads} {input} {params.basename}
        '''

rule hisat2Build:
    input: "genome/genome.fa",
        "genome/"+GTF_FILE
    params:
        basename="indexes/hisat2/reference"
    output:
        output1="indexes/hisat2/reference.1.ht2",
        output2="indexes/hisat2/reference.2.ht2",
        output3="indexes/hisat2/reference.3.ht2",
        output4="indexes/hisat2/reference.4.ht2",
        output5="indexes/hisat2/reference.5.ht2",
        output6="indexes/hisat2/reference.6.ht2",
        output7="indexes/hisat2/reference.7.ht2",
        output8="indexes/hisat2/reference.8.ht2",
    threads: 20
    shell:
        '''
        module load HISAT2
        hisat2_extract_splice_sites.py genome/{GTF_FILE} > genome/genome.ss
        hisat2_extract_exons.py genome/{GTF_FILE} > genome/genome.exon
        hisat2-build -f -p {threads} {input} --ss genome/genome.ss --exon genome/genome.exon {params.basename}
        rm genome/genome.ss
        rm genome/genome.exon
        '''

rule clean:
    shell:
        '''
        rm genome/*
        rm indexes/bowtie2/*
        rm indexes/hisat2/*
        '''

And this is my driver script which I'm submitting as a job itself:
#!/bin/bash

module load snakemake

snakemake --jobs 1 --local-cores 20 --cluster-config /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/cluster.json --latency-wait 120 --max-jobs-per-second 1 -c "sbatch -A {cluster.account} -p {cluster.partition} -n {-cluster.n} -t {cluster.time} -C {cluster.C}" -s indexMake

But I'm getting this error:
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 1
Job counts:
    count   jobs
    1   all
    1   bowtie2Build
    1   getGenome
    1   hisat2Build
    4

rule getGenome:
    output: genome/genome.fa, genome/Mus_musculus.GRCm38.94.gtf
    jobid: 2

RuleException in line 30 of /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/indexMake:
NameError: The name '-cluster' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}

}

I really am not sure what I'm doing wrong here ? Is there a conflict with variables or there is something wrong with the snakemake version I'm using (4.5.0).

(vote to leave open) I know that Snakemake can be used for non-bioinformatics problems. But in practice it is used solely by bioinfo community therefore I think debugging snakefiles is appropriate for this Q&A site.

(vote to close) The official Snakemake docs point folks to Stack Overflow. I suggest we use StackExchange's migrate mechanism to send generic Snakemake questions to SO to keep it all in one place.

I'm voting to close this question as off-topic because it belongs on StackOverflow (although I lack the reputation to initiate the migration).

@DanielStandage it does seem to be aimed primarily at bioinformaticians, and is indeed a tool commonly used by the community, so it seems reasonable that we should support it.

In your submission script you have -n {-cluster.n}, hence the error from snakemake I guess. Are you sure it is not a typo and it should be -n {cluster.n} instead?

Yes, that was indeed the error. I'm not sure how it crept in - working on Sundays is not that productive -


text:
blasting a refseq protein does not show the protein in the result set
Can anyone explain me, why I don't find a specific protein with a blast that was took before from the NCBI refseq database?
Specifically, I was trying to blast the protein with the accession number "NP_420767" and its sequence, respectively, however that protein does not show up in the results. It not only happens when the standard options are chosen, but also, when "Reference proteins (refseq_protein)" as database in the blast options is selected.
I am really puzzled by that.. Shouldn't blast show the initial refseq entry in the results list as it is part of the refseq database and has the same sequence?

Relevant: Duplicate long hits from PSI-BLAST

NP_420767 is represented by the non-redundant refSeq protein WP_010919826, which has the same amino acid sequence. This is not very clearly annotated, but if you scroll down to the sequence in the GenPept entry for NP_420767, you'll see the following:
CONTIG      join(WP_010919826.1:1..799)

See here for more info.

Ah, I understand.. And in the list after clicking on "identical proteins" on the "WP_010919826.1" site, "NP_420767" gets listed.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Getting protein FASTA sequence based on keyword with python
I would like to gather proteins FASTA sequence from Entrez with python 2.7. I am looking for any proteins that have the keywords: "terminase" and "large" in their name. So far I got this code:
from Bio import Entrez
Entrez.email = "example@example.org"

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)
searchResult = Entrez.read(searchResultHandle)
ids = searchResult["IdList"]

handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
record = handle.read()

out_handle = open('myfasta.fasta', 'w')
out_handle.write(record.rstrip('\n'))

However, it can get me several terminases from various organisms, while I need only terminase form bacteriophages (specificly Viruses [taxid 10239], host: bacteria). I've managed to get the nuccore accession ids from NCBI of the viruses I am intersted in, but I don't know how to combine those two informations. The id file looks like this:
NC_001341
NC_001447
NC_028834
NC_023556
...

Do I need to access every gb file of every ID and search for my desired protein in it?

Welcome to Bioinformatics! Did you tried to add the " [taxid 10239], host: bacteria" in the search in Entrez? From the retrieved fastas you could also check if the id is of interest for you and then omit those who aren't

I tried, but I don't know if my syntax is correct. Should it be: term="terminase large bacteria[organism]"

You the term search would be the same query you set on the NCBI page, that's '("Viruses"[Organism]) AND (host[All Fields] AND ("Bacteria"[Organism] OR "Bacteria Latreille et al. 1825"[Organism] OR bacteria[All Fields]))'

Found what I was looking for. In:
searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)

I've added:
searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)

which norrowed down my searches to the desired viruses. Granted it's not filtered by host, but by a taxonomy group, but it is enough for my work.
Thank you @Llopis for additional help

Please remember to accept this answer. The system won't let you accept your own answer just after posting it, but you should be able to now.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Difference between samtools mark duplicates and samtools remove duplicates?
What is difference between samtools mark duplicates and remove duplicates ?  Is it necessary to mark duplicates before removing duplicates with samtools?

samtools rmdup and samtools markdup -r do the same thing. Without the -r flag samtools markdup only flags the duplicates.
You'll have to run samtools fixmate -m and sort the output to add ms and MC tags to prepare the file for markdup..
According to the current documentation, rmdup is obsolete. Please do not use rmdup.

TL;DR: just use markdup.
rmdup removes duplicates from BAM, while markdup, like Picard's MarkDuplicates, marks duplicates by default without hard removal – the latter is usually the desired behavior. In addition, markdup implements a better algorithm that takes care of more corner cases and gives more consistent results.
PS: to emphasize – I see no point to use samtools rmdup nowadays. It has been declared obsolete as @kohlkopf said.

Could you elaborate regarding why not performing hard removal "is usually the desired behavior"? For instance, what useful stuff is typically done using the duplicates?


text:
Issues performing variant calling with GATK
I am trying to perform variant calling on a BAM file generated through STAR version STAR_2.6.0b for wheat genome using GATK haplotypecaller as follows:
gatk HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
02:34:13.450 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
02:34:13.600 INFO  HaplotypeCaller - ------------------------------------------------------------
02:34:13.600 INFO  HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.0.4.0
02:34:13.600 INFO  HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
02:34:13.601 INFO  HaplotypeCaller - Initializing engine
02:34:40.721 INFO  HaplotypeCaller - Done initializing engine
02:34:48.175 INFO  HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
02:34:48.598 INFO  HaplotypeCaller - Shutting down engine
[June 7, 2018 2:34:48 AM PKT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.59 minutes.
Runtime.totalMemory()=10305929216
java.lang.IllegalArgumentException: samples cannot be empty

When I searched for this error the most proposed solution was to validate BAM file as suggested here. The ValidateSamFile gave following output for the bam file:
gatk ValidateSamFile -I sorted.bam -M SUMMARY
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar ValidateSamFile -I sorted.bam -M SUMMARY
02:53:27.498 INFO  NativeLibraryLoader - Loading libgkl_compression.so   from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 02:53:27 PKT 2018] ValidateSamFile  --INPUT sorted.bam --MODE   SUMMARY  --MAX_OUTPUT 100 --IGNORE_WARNINGS false --VALIDATE_INDEX true --INDEX_VALIDATION_STRINGENCY EXHAUSTIVE --IS_BISULFITE_SEQUENCED false --MAX_OPEN_TEMP_FILES 8000 --SKIP_MATE_VALIDATION false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 02:54:08     SamFileValidator        Validated Read    10,000,000 records.  Elapsed time: 00:00:36s.  Time for last 10,000,000:   36s.  Last read position: TGACv1_scaffold_080127_1DS:28,522
INFO    2018-06-07 02:54:44     SamFileValidator        Validated Read    20,000,000 records.  Elapsed time: 00:01:12s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_602746_7DL:127,507
INFO    2018-06-07 02:55:20     SamFileValidator        Validated Read    30,000,000 records.  Elapsed time: 00:01:48s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_641946_U:17,161
INFO    2018-06-07 02:55:55     SamFileValidator        Validated Read    40,000,000 records.  Elapsed time: 00:02:23s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_471755_6AL:43,097
INFO    2018-06-07 02:56:31     SamFileValidator        Validated Read    50,000,000 records.  Elapsed time: 00:02:59s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_434075_5DL:29,039
INFO    2018-06-07 02:57:07     SamFileValidator        Validated Read    60,000,000 records.  Elapsed time: 00:03:34s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_062491_1DL:40,466
INFO    2018-06-07 02:57:42     SamFileValidator        Validated Read    70,000,000 records.  Elapsed time: 00:04:10s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_558158_7AL:13,072
INFO    2018-06-07 02:58:18     SamFileValidator        Validated Read    80,000,000 records.  Elapsed time: 00:04:46s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_435811_5DL:29,008
INFO    2018-06-07 02:58:53     SamFileValidator        Validated Read    90,000,000 records.  Elapsed time: 00:05:21s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_572030_7AS:9,114

## HISTOGRAM    java.lang.String
Error Type      Count
ERROR:MATE_NOT_FOUND    107771
ERROR:MISSING_READ_GROUP        1
WARNING:MISSING_TAG_NM  96856336
WARNING:RECORD_MISSING_READ_GROUP       96856336

In order to remove the errors shown by ValidateSamFile I used FixMateInformation as suggested here. as shown below(sharing part of file):
gatk FixMateInformation -I sorted.bam -O mate.bam -MC
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar FixMateInformation -I sorted.bam -O mate.bam -MC
06:13:13.934 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 06:13:19 PKT 2018] FixMateInformation  --INPUT sorted.bam --OUTPUT mate.bam --ADD_MATE_CIGAR true  --ASSUME_SORTED false --IGNORE_MISSING_MATES true --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 06:13:24     FixMateInformation      Sorting input into queryname order.
INFO    2018-06-07 06:23:26     FixMateInformation      Sorting by queryname complete.
INFO    2018-06-07 06:23:26     FixMateInformation      Output will be sorted by coordinate
INFO    2018-06-07 06:23:26     FixMateInformation      Traversing query name sorted records and fixing up mate pair information.
INFO    2018-06-07 06:23:33     FixMateInformation      Processed      1,000,000 records.  Elapsed time: 00:00:06s.  Time for last 1,000,000:    6s.    Last read position: TGACv1_scaffold_361122_4DS:70,952
INFO    2018-06-07 06:23:42     FixMateInformation      Processed     2,000,000 records.  Elapsed time: 00:00:15s.  Time for last 1,000,000:    8s.    Last read position: TGACv1_scaffold_044571_1BL:711
INFO    2018-06-07 06:23:50     FixMateInformation      Processed     3,000,000 records.  Elapsed time: 00:00:23s.  Time for last 1,000,000:    8s.   Last read position: TGACv1_scaffold_433135_5DL:36,647
 .
 .
 .
 .
 [Thu Jun 07 06:40:50 PKT 2018] picard.sam.FixMateInformation done.    Elapsed time: 27.62 minutes.
 Runtime.totalMemory()=8004304896
 Tool returned:
 0

After performing FixMateInformation I again validated the BAM file using ValidateSamFile but the errors are still the same as shown above which shows Mate not found, which means error is not fixed.
How these errors can be fixed? So that I can perform variant calling using BAM file?

The primary error isn't the missing mate, but likely the missing read groups. Use Picard AddOrReplaceReadGroups to ensure your BAM files have read groups and then run those through the haplotype caller.

Should I run AddOrReplaceReadGroups with default parameters as discussed here or there is something that needs to be extracted and provided as parameters for AddOrReplaceReadGroups?

There are no defaults for the read group information. You can just make up values to put there (usually I just use the sample name). Note that you only need the required fields.

You could add read groups:
['STAR', '--twopassMode', 'Basic', '--genomeDir', '/data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/',
         '--runThreadN', '10', '--limitSjdbInsertNsj', '5000000', '--outSAMtype', 'BAM', 'SortedByCoordinate', '--twopass1readsN', '-1',
         '--sjdbOverhang', '100', '--readFilesIn', sra1 + '_1.fastq', sra1 + '_2.fastq',
         '--outSAMattrRGline', 'ID:' + rg1, 'SM:' + rg1, 'PL:ILLUMINA'], shell=False)

In other words, "add read groups".

@user3018491 on which version of STAR you have used the parameter SortedByCoordinate? As discussed here I tried to use this on version 2.6.0b but it failed.

STAR
Usage: STAR  [options]... --genomeDir REFERENCE   --readFilesIn R1.fq R2.fq
Spliced Transcripts Alignment to a Reference (c) Alexander Dobin, 2009-2015
versions
versionSTAR             020201

The problem is that you are missing read groups. You can add them directly to bam files you already have using picard tools. Something like
java -jar ./software/picard.jar AddOrReplaceReadGroups \
  I=input.bam \
  O=output.bam \
  RGLB=lib1 \
  RGPL=illumina \
  RGPU=unit1 \
  RGSM=hello

Hello 罗素幸福之路, it would be nice to explain what the code you have posted means. I tried to add at least one sentence, feel free to [edit] it and improve it further.


text:
Disk space error while aligning reads using STAR
Hi am trying to align RNA-seq reads using STAR through following commands,
STAR --runMode alignReads --runThreadN 26  --outSAMtype BAM  SortedByCoordinate Unsorted --sjdbOverhang 100  --genomeDir wheat --readFilesIn G3_cleaned_R1.fastq G3_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3
Jun 05 13:53:59 ..... started STAR run
Jun 05 13:53:59 ..... loading genome
Jun 05 14:07:12 ..... started mapping
Jun 05 18:42:47 ..... started sorting BAM

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin   does not agree with the actual size on disk: 509186583   0   45

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 648593301   0   47

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 660412146   0   46
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 470931698   0   44

Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 653155717   0   43

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Segmentation fault (core dumped)

Although there is much free space on my system as shown below:
> df -Th
  Filesystem     Type      Size  Used Avail Use% Mounted on
  devtmpfs       devtmpfs   63G     0   63G   0% /dev
  tmpfs          tmpfs      63G  136K   63G   1% /dev/shm
  tmpfs          tmpfs      63G  1.9M   63G   1% /run
  tmpfs          tmpfs      63G     0   63G   0% /sys/fs/cgroup
  /dev/sda3      ext4      3.6T  1.4T  2.2T  39% /
  /dev/sda2      ext4      9.8G   57M  9.2G   1% /boot
  tmpfs          tmpfs      13G     0   13G   0% /run/user/480
  tmpfs          tmpfs      13G   16K   13G   1% /run/user/1000

How this problem can be solved?

Are you an older version of the STAR aligner? A similar error has been discussed here: https://groups.google.com/forum/#!topic/rna-star/OTUG26gKKmo

You are using 26 threads, while you have 125GB of RAM, this means per thread you have 4.8 GB RAM. Your genome size is 13 GB (according to your previous post). An answer to your previous post suggested to reduce the number of threads. Did you try to run with lower number of threads?

@b.nota I ran the job discussed in previous post with same number of threads.

@Kasper I am using latest version of STAR

Report the bug on github and sort with samtools.

@DevonRyan ok I will post it there. Can I use --limitBAMsortRAM to solve this problem?

@AmmarSabirCheema I don't think so, since any time you need to increase that there's an error message saying so. I think this is just a bug somewhere. Just sort the unsorted output with samtools or sambamba.

STAR versions 2.6.0b and 2.6.0a are unstable,  as written Issues performing variant calling with GATK you are using version 2.6.0b. You should switch to version 2.6.0c which is stable.  Also avoid to perform multiple alignments in the same directory.


text:
Reject reads with low quality bases from a Bam file through pysam
I have a code below:
def FilterReads(in_file, out_file):

def read_ok(read):
    """
    read_ok - reject reads with a low quality (<5) base call
    read - a PySam AlignedRead object
    returns: True if the read is ok
    """
    if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
        return False
    else:
        return True

_BASE_QUAL_CUTOFF = 30

bam_in = pysam.Samfile(in_file, 'rb')
bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)

for read in bam_in.fetch():
    if read_ok(read):
        bam_out.write(read)

This code works fine by

Rejecting reads having a base with phred quality score below five
But it first takes a BAM file
And then creates a filtered BAM file for further analysis but this takes a lot of time.

So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it?
Or can I modify this code to perform the same function

Why are you sorting the file prior to iterating over it?  Your function doesn't require it to be sorted and since you are not using the random indexing pysam doesn't require it either.

if …: return False else: return True is an anti-pattern. Just write return not …. Or, better yet, logically invert the condition; in your case, you could write return all([ord(c) - 33 >= _BASE_QUAL_CUTOFF for c in list(read.qual)]) (this follows from straightforward application of de Morgan’s laws).

This is probably not the cause of the slowdown, but you can avoid two list creations per read by refactoring line 9 as any(ord(c)-33 < _BASE_QUAL_CUTOFF for c in read.qual):

I don't think pileup is going to save you any time here, as it has to load the same data that you are loading with the fetch method.

For can do this by accessing the basecall qualities from PileupRead.alignment. For example:
_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

This involves checking every read at every position. You could memoise this task:
def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

Or, using functools.lru_cache (thanks @winni2k):
from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

FYI: There's a decorator in the library that does this and more! https://docs.python.org/3/library/functools.html#functools.lru_cache


text:
Increase number of threads for GATK 4.0 HaplotypeCaller
I am using GATK version 4.0, I tried to use multiple threads for calling variants using HaplotypeCaller using following command
gatk --java-options -Xmx90G  -nt 28 HaplotypeCaller -I output.bam -R wheat_ref.fa -O final.vcf

and the error is
'-nt' is not a valid command

I have also used --num_threads but again the error was the same, as written here this parameter is available in previous versions of GATK but I don't know what's its status in current version of GATK i.e GATK 4.0. Whether it is deprecated or replaced by some other parameter?
Is there any way of calling multiple threads for variant calling with GATK 4.0?

I am assuming you are using the non-SPARK version of the method. To specify the number of threads you wish to use with HaplotypeCaller, include --native-pair-hmm-threads (documentation). This will only parallelize the pair hidden Markov models (pair HMM) process. This is an argument for the method HaplotypeCaller itself, so it should be located after the method is called as opposed to one of the --java-options.
I do recommend trying out the SPARK implementations of the GATK tools. The whole method will take advantage of the threads you specify. I've been running SPARK locally (i.e. not on Firecloud) on computational servers. I've been impressed with the enhanced performance.

I used this command now gatk HaplotypeCaller --native-pair-hmm-threads 28  -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf but it took same time and also I specified 28 threads and 28 was shown in requested threads but system monitor was showing only 3 ~ 4 threads are used as shown here https://imgur.com/a/JJJYmB7.


text:
Filtering bases based on phred qualities with pysam
Is there a way to filter bases in BAM files based on phred quallities through python's pysam ?
I have a code here that

Takes the nucleobases per position from a BAM file using pysam's pileup function
Saves it in ReverseList and ForwardList based on both strands (i.e forward and reverse), I want to reject those bases that have phred quality below 25 to be not stored in the ForwardList and ReverseList lists so that they are not used for further analysis.
 samfile = pysam.Samfile( filename, "rb" )
 ReverseList = [''] * lenref
 ForwardList = [''] * lenref
 for pileupcolumn in samfile.pileup() :
     for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
               continue
         if not pileupread.is_del and not pileupread.is_refskip:
               if pileupread.alignment.is_reverse: #negative
                      ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
               else:
                      ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

 samfile.close()

Where lenref = 16569(length of mitochondrial genome) and filename is the name of BAM file. I want to filter based of phred qualities of bases.

Am I the only one who is missing the point of filtering alignments by phred scores? Phred scores are assigned during base calling, if you want to filter them out, why after mapping but not before? Also variant callers take into account phred scores, having low phred scores should not worsen quality of calling variants.

I am making a tool for analyzing BAM file so I want to filter it again to improve my results.

It is also reported in previous studies.

Welcome to Bioinformatics! This question seem to need more information to get a good answer (Although the existing one answers how to filter as required). Why do you want to filter the bases? What are you trying to achieve?

The PileupRead object as a query_position attribute, which you can use for this:
for pileupcolumn in samfile.pileup() :
    for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
             continue
         if not pileupread.is_del and not pileupread.is_refskip:
             if pileupread.alignment.query_qualities[pileupread.query_position] < 10:
                 # Skip entries with base phred scores < 10
                 continue
             if pileupread.alignment.is_reverse: #negative
                   ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
             else:
                   ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

Note the 6th line, which implements a filter with a threshold of 10.


text:
snakefile not able to find files for rule trimming
I am trying to run a snakemake pipeline for small-RNA samples and having  problems. The snakefile is given below,
configfile: "config.yaml"
#include: "utils.py"
#Command line syntax:
#snakemake -s snakefile.smk --cores 8 --config prefix= lib=
class MyException(Exception):
    def __init__(self, message):
        self.message = message
        super().__init__(self.message)

with open(config['IDs'], 'r') as file:
    # Read all the lines of the file into a list
    lines = file.read().splitlines()
rule all:
    input:
        expand("fastq/{prefix}.fastq.gz", prefix=lines),
        expand("trimfastq/{prefix}/{prefix}_trimmed.fq.gz", prefix = lines),
        expand("trimfastq/{prefix}", prefix = lines),
        expand("reportfalco/{prefix}", prefix=lines),
        expand("sam/{prefix}Aligned.out.bam", prefix = lines),
        expand("sam/{prefix}.s.bam", prefix=lines),
        expand("sam/{prefix}.s.bam.bai", prefix=lines)

#Checking fastq quality
rule testquatility:
    input:
        fq = "fastq/{prefix}.fastq.gz"
    output:
        falrep = directory("reportfalco/{prefix}")
    shell:
        """
        /home/tools/Falco/falco/bin/falco {input.fq} -o {output.falrep}
        """
#Trimming
rule trimming:
    input:
        fq_in = "fastq/{prefix}.fastq.gz"
    output:
        trimdir = directory("trimfastq/{prefix}"),
        fq_out = "trimfastq/{prefix}/{prefix}_trimmed.fq.gz"
    run:
        if config['Lib']=='nextflex':
            shell("/home/tools/TrimGalore-0.6.6/trim_galore -a TGGAATTCTCGGGTGCCAAGG --clip_R1 4 --three_prime_clip_R1 4 --length 18 {input.fq_in} --output_dir {output.trimdir}")
        elif config['Lib']=='comboseq':
            shell("/home/tools/TrimGalore-0.6.6/trim_galore -a AAAAAAAA --length 15 {input.fq_in} --output_dir {output.trimdir}")
        elif config['Lib']=='trilink':
            shell("/home/tools/TrimGalore-0.6.6/trim_galore -a TCACCGACTGCCATAGAG --small_rna {input.fq_in} --output_dir {output.trimdir}")
        else:
            raise MyException("Library strategy is expected to be either: trilink, comboseq or nextflex; but is" + str(config['Lib']))
            with MyException as e:
                print("\033[91mError when we try to execute trimmimg :", e.message, "\033[0m")

#Alignment
rule align:
    input:
        ref="/references/DERFINDER_ref/STAR_INDEX",
        trim_in="trimfastq/{prefix}/{prefix}_trimmed.fq.gz"
    params:
        pr="sam/{prefix}"
    output:
        "sam/{prefix}Aligned.out.bam"
    run:
        if config['Lib']=='nextflex' or config['Lib']=='trilink':
            shell("/home/tools/STAR-2.7.9a/bin/Linux_x86_64/STAR --runThreadN 16 --genomeDir {input.ref} --readFilesIn {input.trim_in} --readFilesCommand gunzip -c --outFileNamePrefix {params.pr} --outFilterScoreMinOverLread 0 --outFilterMatchNmin 16 --outFilterMatchNminOverLread 0 --outFilterMismatchNoverLmax 0.025 --alignIntronMax 1 --alignEndsType EndToEnd --outSAMtype BAM Unsorted")
        elif config['Lib']=='comboseq':
            shell("/home/tools/STAR-2.7.9a/bin/Linux_x86_64/STAR --runThreadN 16 --genomeDir {input.ref} --readFilesIn {input.trim_in} --readFilesCommand gunzip -c --outFileNamePrefix {params.pr} --outFilterScoreMinOverLread 0 --outFilterMatchNmin 15 --outFilterMatchNminOverLread 0.9 --outFilterMismatchNmax 1 --outFilterMismatchNoverLmax 0.3 --alignIntronMax 1 --outSAMtype BAM Unsorted")
        else:
            raise MyException("Library strategy is expected to be either: trilink, comboseq or nextflex; but is" + str(config['Lib']))
            with MyException as e:
                print("\033[91mError when we try to execute trimmimg :", e.message, "\033[0m")

#Housekeeping
rule sort:
    input:
        al_in="sam/{prefix}Aligned.out.bam",
    output:
        sambam="sam/{prefix}.s.bam"
    shell:
        """
        samtools sort -@ 16 -o {output.sambam} {input.al_in}; samtools index {output.sambam}; rm {input.al_in}
        """

#index reads
rule index:
    input:
        al_in="sam/{prefix}.s.bam",
    output:
        sambam="sam/{prefix}.s.bam.bai"
    shell:
        """
        samtools index {input.al_in}
        """

The config.yaml file is given below,
  IDs: IDs.txt
  Lib: 'comboseq'

The file IDs.txt is given below,
  sample.fastq.gz

When I run the above snakemake pipeline using conda environment and the command given below,
  snakemake --cores 8 --set-threads align=16 sort=16 index=16

then it gives the error given below,
  MissingInputException in rule trimming in file /mnt/data/Snakefile,
  line 34:
  Missing input files for rule trimming:
  output: trimfastq/sample.fastq.gz,
  trimfastq/sample.fastq.gz/sample.fastq.gz_trimmed.fq.gz
  wildcards: prefix=sample.fastq.gz
  affected files:
  fastq/sample.fastq.gz.fastq.gz

Do you have any solution for this problem?
Thanks,

The Problem is, that your ID is called sample.fastq.gz and you tell snakemake to look for {ID}.fastq.gz. Therefor you get the error sample.fastq.gz.fastq.gz can not be found.
Thus you have two options. Either remove the .fastq.gz-suffix in your IDs.txt or remove it in the snakefile.


text:
Is there a standard k-mer count file format?
I am doing a research project involving calculating k-mer frequencies and I am wondering if there is any standard file format for storing k-mer counts.

Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the dump command (but its default output format is a custom binary format):

The dump subcommand outputs a list of all the k-mers in the file associated with their count. By default, the output is in FASTA format, where the header line contains the count of the k-mer and the sequence part is the sequence of the k-mer. This format has the advantage that the output contains the sequence of k-mers and can be directly fed into another program expecting the very common FASTA format. A more convenient column format (for human beings) is selected with the -c  switch.

Jellyfish changed their internal format between v1 and v2 (both not FASTA), because they changed to doing counts based on bloom filters. Jellyfish2 has an optional two-pass method that sets up a bloom filter intermediate file to record kmers, and multiple different final reporting formats.
Khmer also uses bloom filters, but in a slightly different way. It also has been extended to be useful for partitioning and comparing datasets.

Based on my experience reviewing k-mer counting software and as a core contributor to the khmer project, I can confidently say that there is no widely used standard format.


text:
Fasta Sequence Identifier format?
I have been analyzing some virus DNA from the NCBI databases. The fasta sequences that I recieve from them have header lines that look like this:
>gi|61393989|gb|AY848686.1|

I think that the second number is the GenBank ID, but I have several questions:
1) What does the "gi" in the first ID stand for?
2) Is this a standard format for storing a series of IDs associated with sequences? That is, should I expect to see this format elsewhere and be prepared to parse fasta headers like this in code that I write?
3) It seems that AY848686 is also a valid GenBank ID, and I see the .1 missing/present in GenBank IDs in various places. What does the .1 mean at the end of GenBank IDs and why is it there/does it need to be there?
Thanks in advance for you answers.

The gi is an abbreviation for "Genbank identifier".
This is a pretty standard convention used by data stored in NCBI databases. There used to be a pretty comprehensive description of the conventions used at NCBI (I wouldn't say it was a standard or specification, just convention) here, but this page is no longer available it seems.
The AY848686 label is an accession value, and the appended number is a version number.

See this page for more info.
The short answer is: these are all internal conventions. None of it needs to be there, except perhaps to ensure NCBI's internal tools work as expected.

UPDATE: I found a non-official copy of the old docs here. Look for the section starting with "FASTA Defline Format". I've copied it here to prevent further problems with bitrot in the future.
FASTA Defline Format
--------------------
The syntax of FASTA Deflines used by the NCBI BLAST server depends on
the database from which each sequence was obtained.  The table below lists
the identifiers for the databases from which the sequences were derived.

  Database Name                         Identifier Syntax

  GenBank                               gb|accession|locus
  EMBL Data Library                     emb|accession|locus
  DDBJ, DNA Database of Japan           dbj|accession|locus
  NBRF PIR                              pir||entry
  Protein Research Foundation           prf||name
  SWISS-PROT                            sp|accession|entry name
  Brookhaven Protein Data Bank          pdb|entry|chain
  Patents                               pat|country|number
  GenInfo Backbone Id                   bbs|number
  General database identifier           gnl|database|identifier
  NCBI Reference Sequence               ref|accession|locus
  Local Sequence identifier             lcl|identifier

For example, an identifier might be "gb|M73307|AGMA13GT", where the "gb" tag
indicates that the identifier refers to a GenBank sequence, "M73307" is its
GenBank ACCESSION, and "AGMA13GT" is the GenBank LOCUS.

"gi" identifiers are being assigned by NCBI for all sequences contained
within NCBI's sequence databases.  The 'gi' identifier provides a uniform
and stable naming convention whereby a specific sequence is assigned
its unique gi identifier.  If a nucleotide or protein sequence changes,
however, a new gi identifier is assigned, even if the accession number
of the record remains unchanged. Thus gi identifiers provide a mechanism
for identifying the exact sequence that was used or retrieved in a
given search.

For searches of the nr protein database where the sequences are derived
from conceptual translations of sequences from the nucleotide databases
the following syntax is used:

                     gi|gi_identifier

An example would be:

        gi|451623           (U04987) env [Simian immunodeficiency..."

where '451623' is the gi identifier and the 'U04987' is the accession
number of the nucleotide sequence from which it was derived.

Users are encouraged to use the '-I' option for Blast output which will
produce a header line with the gi identifier concatenated with the database
identifier of the database from which it was derived, for example, from a
nucleotide database:

        gi|176485|gb|M73307|AGMA13GT

And similarly for protein databases:

        gi|129295|sp|P01013|OVAX_CHICK

The gnl ('general') identifier allows databases not on the above list to
be identified with the same syntax.  An example here is the PID identifier:

        gnl|PID|e1632

PID stands for Protein-ID; the "e" (in e1632) indicates that this ID
was issued by EMBL.  As mentioned above, use of the "-I" option
produces the NCBI gi (in addition to the PID) which users can also
retrieve on.

The bar ("|") separates different fields as listed in the above table.
In some cases a field is left empty, even though the original
specification called for including this field.  To make
these identifiers backwards-compatiable for older parsers
the empty field is denoted by an additional bar ("||").

Thanks for the quick and accurate answer. As a follow up question- is there some reason why NCBI has accession numbers as well as an ID numbers?

@JonDeaton history, mostly. The GIs are far older. Also see the "Historical Note" section in the link Daniel gave you.

@JonDeaton For parsing be aware that NCBI have phased out GI numbers from fasta headers since Sept 2016. I couldnt find a new source for this info but please see this note from their blog https://ncbiinsights.ncbi.nlm.nih.gov/2016/07/15/ncbi-is-phasing-out-sequence-gis-heres-what-you-need-to-know/

That is a GenBank gene ID and it (apparently) stands for "GenInfo Identifier". The full GenBank ID is gi|61393989.

Yes, this is the standard format for sequences provided by NCBI. However, you can have all sorts of extra details appended to that. They (NCBI) generally use | as a field separator, so you could try to write code that parse that, but I really wouldn't recommend it. The FASTA headers have traditionally been free. People can put anything they like there apart from a \n, so I really wouldn't put any effort into writing code to parse these headers.
If you need a short name for the sequence, just take the first N characters or, if it matches ^..|.+?|, take the characters until the second |. Seriously though, don't bother parsing this it is very rarely helpful and is almost certain to break as soon as you use a different database.

The AY848686.1 is a GenBank accession. The gb indicates that it's a GenBank accession and the .1 means it is the 1st version of that accession.

As for why there are two (source):

Unlike the gi number system, in which sequence identification numbers were not necessarily consistent across the databases (e.g., GenBank and EMBL could each assign their own gi number to a sequence), the new system is designed to ensure consistency. It is also designed to show a relationship between a sequence identification number and the accession number of the record in which it is found. In contrast, GI numbers are assigned consecutively and bear no resemblance to the accession number. Finally, the new system allows the assignment of alphanumeric protein IDs to proteins translations within nucleotide sequence records. The protein IDs contain three letters followed by five digits, a period, and a version number.


text:
How to convert files to ADAM format?
I would like to convert BAM and VCF files to ADAM format.
How do I do that?

These transformations are quite simple with the adam-submit script packaged with ADAM.
Transformation of BAM files
adam-submit -- transformAlignments sample.bam sample.alignments.adam

Transformation of VCF files
adam-submit -- transformVariants sample.vcf.gz sample.variants.adam

Could you provide a link to the package? It seems like it is this repository but I'm not sure.


text:
Bismark report before and after Cytosine Methylation Extraction
Looking at the sample report from Bismark webpage, there are two sections: Cytosine Methylation and Cytosine Methylation after Extraction. However even after reading User Guide, I still cannot understand the difference between these two sections.
As far as I understand this, Cytosine Methylation section shows summary statistics about Cytosine Methylation in sequenced data and Cytosine Methylation after Exctraction shows summary statistics after a Methylation calling, but is not the Methylation calling already part of the first section?
My question therefore is: What is the difference between these two sections of Bismark report?

Did you try contacting the authors of Bismark? It doesn't seem to be highly used, so I doubt we have here someone who could answer you. (If you find the answer, then post it here for others please)

I was also looking for an answer to the same question and found this issue. have a look at this post https://github.com/FelixKrueger/Bismark/issues/245


text:
Standard Cutoff for Moderated T-statistics
I'm looking at some microarray data. For the first time I've calculated a moderated T statistic from limma.
Is there any standard practice for where to cut off that value? For log2 fold change I usually cut off at +/- 1.5 and adjust accordingly; for adjusted p value I cut off at p<0.01 and adjust accordingly.
These moderated T statistics seem a little higher than log2 fold change and I'm not sure what is generally accepted as significant.
Thanks!

If you're filtering by the (adjusted) p-value, why do you want to also filter by moderated T-statistic? The p-value basically comes from the moderated T-statistic, after all.

I would like to have values that are real but are also above a particular magnitude. The P value gives me an idea of whether the finding is real, but I also need a way to cut-off on magnitude. I'd usually do this by filtering l|og2Fold|>1.5, but I was hoping to have a similar moderated T-statistic cutoff. Or am I misinterpreting the moderated T statistic?

You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (adjusted) p-value and fold-change thresholds.


text:
What are the count units in DEseq2?
This is a pretty silly/simple question, I suppose.
What units are DESeq2 counts? Or are the units arbitrary but internally consistent estimate of actual reads?

The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.


text:
RNAseq: Z score, Intensity, and Resources
I'm very new to bioinformatics in general, and I'm trying to understand some basic concepts.
I have RNAseq data, and bioinformatics people tell me that intensities cannot be compared across patients. So there are all of these pipelines to compare intensities to Z scores--are those as simple as just plugging data into a bioconductor package?
It would be great to get some overview description of why Z values are important/why you can't compare intensities across patients, and/or some pointers towards any resources that I can read to learn more.
Thanks much for your time!

Your wording suggests that you have microarray data rather than RNAseq, is this correct? Z values are rarely used in RNAseq, rather counts per gene/feature or TPMs are typically used.

I have some of both--even though I think I've seen Z-scores from TCGA illumna sets? Or perhaps I was mistaken.

It depends on what test or analysis you want to do, whether you need intensities (expression values) or z-scores.
If you want to do statistical analysis, such as finding differentially expressed genes between groups of patients (e.g., with limma), you don't want to use z-scores. But you use normalized intensities (for microarrays) or start with raw counts for RNAseq. Then follow the user guide (limma has a great user guide, also for starters).
For visualization in heatmaps or for other clustering (e.g., k-means, fuzzy) it is useful to use z-scores. Z-scores are a form of transformation (scaling), where every genes is sort of "reset" to the mean of all samples, using also the standard deviation. If you want to know exactly what a z-score is, a simple google search can tell you the details.
In R you can use the scale function for z-score transformation. Be aware that the function works on columns though. Which means you have to transpose your matrix first if genes are in the rows, and then transpose them back after transformation.

Rather than transposing the matrix I’d suggest using apply with scale.

In cancer genomics experiments often involve sequencing a large number of tumour samples, but few or no matched normals. This is partly due to financial constraints, but also due to ethical ones: it is easy to get ethical approval for tumour tissue, which is removed as a matter of course in a patient's care anyway, whereas it is hard to get approval for normal tissue, which would not normally be removed and would require an additional procedure.
However, this presents a problem. Most of the procedures we use in standard gene expression analysis involve comparing one group to another after suitable normalisations of the data. But here we have no suitable comparison group. To get around this problem in cancer RNA-seq analysis we do something called outlier detection. That is, we assume that for a given gene, the majority of the patients have "normal" expression, and look for the few patients that don't look like the rest. We do this with the Z-score.
To calculate the Z-score for an observation we subtract the mean of all observations and divide by the standard deviation. Thus the Z score of an observation is how many standard deviations an observation is from the mean of all observations - or how unusual it is.
Thus for Gene A in Patient 1, we calculate how many standard deviations Gene A is from the mean of Gene A across all patients. A very big (or small/negative) value tells us that the expression of Gene A is unusual in patient 1 compared to the other patients.
(NB: in some cancer studies they use the median and median absolute deviation rather than the mean and standard deviation to compute a "robust Z-score" as the data is unlikely to be normally distributed).

Interesting and well explained. Thanks. About the use of the median, I would add that this makes the results less sensitive to outliers.


text:
Chip-Seq data description
I have the following Chip-Seq data and could not found a description of it. Can you help me with it? More, I would like to find out the order of the nucleotides in binding regions; that is, what number they are.
@GWZHISEQ01:319:HA1NPADXX:1:1101:1182:2196 1:N:0:ACTGAT
GGAGGAAGTGTGTTGCTGCAGACAGTAGCAACCAGACCCACACTGCGCGTA
+
CCCFFFFFHHFHHIJJJJIJJIIIIHGJJJJJJJCGHIJJIJJFHIIJJII
@GWZHISEQ01:319:HA1NPADXX:1:1101:1272:2109 1:N:0:ACTGAT
NCCACGCTAGGCTCAGCTTGTCGGCCTGGCTAAGCAGTTGCGAAAGTGCGC
+
#1=DDDFFHHHHHJJIJJJJJJIGEHGGJIGHJFHGJIJJJJJJJ@CHCBG
@GWZHISEQ01:319:HA1NPADXX:1:1101:1418:2141 1:N:0:ACTGAT
GCACGCACTACCCAGAGATCATCCAAAGCCTGAAGCCACAGGGCGCACTCG
+
CCCFFFFFGHGHHJJJIIJJJJIIHGGGIJIIJJIIIJJJIJJJIJIJIJH
@GWZHISEQ01:319:HA1NPADXX:1:1101:1302:2160 1:N:0:ACTGAT
ACACTCTTTCCCTACACGACGCTCTTCCGAGATCGGAAGAGCACACGTCTG
+
CCCFFFFFHDHHHJJJJIJJJJIIIIIJIHJIJJJIIJHIIJJJJJJJJHH

Cross-posted on biostars

Is that a problem?

It's frowned upon at the very least.

I would argue that cross-posting is not a problem provided you reference the cross-posts on each site.

That's a fastq file, you will want to align it to the genome, call peaks, and then use something like MEME to determine binding motifs.

how do I align it to the genome?

@frosh You need more help than a Q&A site can offer. I recommend that you attend a workshop on NGS data analysis, or find a collaborator who can do the analysis for you.

You're right. I will delete the question.

You don't need to delete it.

The following mini-courses may be useful: http://rafalab.github.io/pages/harvardx.html


text:
Identity plus Similarity Interpretation
I am helping a colleague to interpret protein multiple sequence alignment results. My colleague used a third-party proprietary solution (I know...) to calculate the alignment, and one of the resulting tables is a pairwise distance matrix with identity % in the top triangle and identity plus similarity in the bottom triangle.
I am struggling to interpret identity plus similarity, despite the fact I found the metric in several papers.

A guess; Identity: the amino acids that are equal. Similarity: The amino acids of the same group.

I guess you already know what identity means. But for completeness, it is the % of pairwise alignment position which have an identical amino acid residue. For the identity+similarity it is important to know that amino acids can be grouped according to their physicochemical properties e.g. charge, size etc. While there are standard groupings for amino acids it may be that the proprietary algorithm uses different groupings. Either way, % identity+similarity would imply the number of pairwise alignment positions which contain either the same amino acid or amino acids from the same group (however they define groups).


text:
Comparing two files with overlapping regions and get their associated information
I have two file that have overlapping regions. One file has chromosomes, position, CNV data and samples.Other file has chromosomes, positions and ensemble ids of genes.
Here is some lines of two files.
contig  start   end     type  bi-allelic  S1  S2  S3
chr1    12048   177033    DUP FALSE   -1  22  24
chr1    12048   89237     DUP FALSE   12  10  -1
chr1    2985712 3355300   DEL   TRUE  2   2   2

Second file is like this
chr1    901877  911245  gene_id "ENSG00000187583.6"
chr1    22013   41670 gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

I want output like this.
chr1    22013   41670 gene_id "ENSG00000160075.10"  DUP FALSE   12  10  -1
chr1    2985732 3355185 gene_id "ENSG00000142611.12" DEL   TRUE  2   2   2

Is there any way to do this?

You can use bedtools intersect and cut the results:
$ cat file1.bed
chr1    12048   177033  DUP FALSE   -1  22  24
chr1    12048   89237   DUP FALSE   12  10  -1
chr1    2985712 3355300 DEL TRUE    2   2   2

$ cat file2.gtf
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

$ bedtools intersect -a file2.gtf -b file1.bed -wb | cut -f 1-4,8-11
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   12  10
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   -1  22
chr1    2985732 3355185 gene_id "ENSG00000142611.12"    DEL TRUE    2   2


text:
Differences between TM and GDT-TS scores for structure comparison
What are the relative advantages / disadvantages of the TM (template modelling) and GDT (global distance test, total score) scores, used in protein structure prediction competitions, such as CASP? When should I choose one over the other?

Hi. This question seems to need more information to get a good answer as it looks (to me) too general. Could you please [edit] adding what have you found/read about each tool, or any paper you read on the topic to help me answer you?


text:
VerifyBamID freemix score
We all have used freemix score at some stage to check contamination or swaps in our sequencing experiment, but can anyone once for all explain how the score is calculated?

Have you read the paper? As an aside, we haven't all used VerifyBamID, it's mostly only useful in human experiments.

@DevonRyan Yes I have read the paper but could not grasp how the freemix score is estimated? Posted the same question on Biostars about a year ago but did not get any answer!!

Here's my current understanding of the freemix metric output by verifyBamID. The freemix score is a measure of contamination that is calculated without knowing the genotypes of the individuals. In the paper it is equation (2), and the authors refer to it as the "sequence-only" method. However, in addition to the sequencing reads, it requires some knowledge of SNPs and their allele frequency from the population. They demonstrate in the paper that using the SNPs and allele frequency estimates from a different population actually reduces the estimated level of contamination.
As for intuition for what equation (2) is measuring, imagine what would happen to the number of heterozygous sites detected if DNA from two different individuals was sequenced together. The total number of heterozygous sites would increase because it would count all sites where the two individuals are homozygous for opposite alleles or at least one of the two individuals is heterozygous. This would lead to an increased fraction of heterozygous sites over that expected from Hardy-Weinberg Equilibrium (i.e. 2pq).
Here's the explanation of freemix directly from the author, Hyun Min Kang:

The key idea of FREEMIX estimate is to use excessive heterozygosity to estimate the level of contamination. Especially for common SNPs, you will observe higher fraction of heterozygous alleles than 2*p*(1-p), and it turns out that you can quantify the contamination very well if you know the population allele frequency already. If you do not have accurate population allele frequency information, than it would be harder to estimate FREEMIX parameters using verifyBamID.

Additional sources:

Verifying Sample Identities - Implementation
Interpreting output files
verifyBamID Google Group


text:
Where to download baseline/average gene expression level of all human coding genes?
I am looking for the most appropriate dataset for downloading baseline gene expression level across all human coding genes during development. I am aware that EMBL Expression Atlas is one of the resources that provide such information, but I am ideally looking for species averaged (and probably scaled) expression data for all human coding genes during early stages of development!
Any insight is much appreciated!

There is not going to be a public database with this data. Apart from anything else, generating data for early human development is difficult and ethically tricky. Also, most people who care about development would probably think that a whole body average would be meaningless.
The best you could do would be to select tissues you were interested in from something like Expression Atlas and calculate the averages yourself.


text:
Get protein names corresponding to PDB ID
I have a list of about 4000 PDB IDs and would like to get the actual names of the proteins (e.g. lactate dehydrogenase, cytochrome c). I tried the batch header section at the Protein Databank Download page but it refused to accept my PDB IDs in formats (xxxx or xxxx.pdb, individually or space-separated) that worked in an interactive search for the protein structure.
Any suggestions?

Have you tried writing a PyMol script?

Perhaps you could use ftp://ftp.wwpdb.org/pub/pdb/derived_data/index/entries.idx which contains HEADER and COMPOUND fields extracted from all PDB entries.

@Michael — Could you clarify please. A script tTo do what on what? And why in PyMol?

@marcin — Thanks, but the problem with the HEADER and COMPOUND fields is that they are such utter unparsable anarchy. I could write a script to get these fields from all my pdb files, but take for example 102L — I would get as HEADER "HYDROLASE(O-GLYCOSYL)" (useless) and as COMPOUND "HOW AMINO-ACID INSERTIONS ARE ALLOWED IN AN ALPHA-HELIX OF T4 LYSOZYME ENTEROBACTERIA PHAGE T4" where the name of the protein — Lysozyme — is buried. The prehistoric format of pdb files is the problem, but surely someone, somewhere has done something with them.

In which case, going from PDB to Uniprot using SIFTS https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html is the best bet as the Uniprot entries will have all the synonyms you could have for.

@MatteoFerla Thanks. Will check that out.

You can use one of the UniProt Protein APIs.
As you said you have your pdb entries in a text file line by line you can, like this example.txt containing:
1brr
4lzm
2dyi

Using the commandline, you can use a little script like this to download the name, if it is available for the given pdb entry.
while read line;
do
 curl -X GET --header 'Accept:application/json' "https://www.ebi.ac.uk/proteins/api/proteins/pdb:$line" |
 jq -r '.[].protein.recommendedName.fullName.value'  |
 sed "s/^/$line\t/" >> pdb_names.txt;
done < example.txt;

You need to have curl, sed and jq installed on your system.
This gives you following output in pdb_names.txt
1brr    Bacteriorhodopsin
4lzm    Endolysin
2dyi    Ribosome maturation factor RimM

Update
if you want to speed it up, you can run it with parallel.
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" |   jq -r ".[]. .protein.recommendedName.fullName.value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

With the -j option you call how many jobs should run in parallel. The limit of the UniProt API is 200 request per second per user.

Update 7. Nov 2020
To get another info beside the protein name, you need to know how the JSON response from UniProt looks like.
To get also the scientific name, you can run following command:
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" | jq -r ".[] | .protein.recommendedName.fullName.value + \" - \" + .organism.names[0].value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

As result you get this:
1brr    Bacteriorhodopsin - Halobacterium salinarum (strain ATCC 700922 / JCM 11081 / NRC-1)
4lzm    Endolysin - Enterobacteria phage T4
2dyi    Ribosome maturation factor RimM - Thermus thermophilus (strain HB8 / ATCC 27634 / DSM 579)

sed should also be available for you. You can install jq and parallel with brew

Many thanks. Got this running. Only problem is that on OS X the \t tab escape does not work (it is apparently non-standard) and none of the answers here worked for me. However it is easy enough to insert a tab after the four-character PDB ID programmatically.

Check this answer https://stackoverflow.com/a/6392271 for your tab problem.

Thanks again. I know it is bad form to ask supplementary questions, but could you possibly provide a link in your answer to the particular page in the API which describes the need to specify all of the four "dotted" terms. I went to https://www.ebi.ac.uk/proteins/api to try to find out how to construct a query to get the species (so I didn't have to ask) but must be looking in the wrong place. I can find an example XML for a protein under proteins/proteins that includes some of these terms, and various unpromising json thingies. NCBI provides intelligible APIs (or used to) — whereas EBI…

I've updated my answer. The dotted format is not a problem of the API. It is just how jq is getting the keys from the JSON formatted result. https://stedolan.github.io/jq/manual/#Basicfilters gives you a basic overview of how to use jq.

Assuming you can use R, have you tried with biomaRt? For example, using 2bhl (my PhD lover :D)
library(biomaRt)
ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl")

# get list of all available info
filters <- listFilters(ensembl)
attributes <- listAttributes(ensembl)

getBM(attributes=c('hgnc_symbol','ensembl_gene_id','entrezgene_id',
                'protein_id','description',"superfamily"),
      filters = 'pdb',
      values = "2bhl",
      mart = ensembl)

Returns
hgnc_symbol ensembl_gene_id entrezgene_id protein_id description
1         G6PD ENSG00000160211          2539   ADO22353 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
2         G6PD ENSG00000160211          2539   CAA27309 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
3         G6PD ENSG00000160211          2539   AAA63175 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
4         G6PD ENSG00000160211          2539   AAA52500 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
5         G6PD ENSG00000160211          2539   AAA52501 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
6         G6PD ENSG00000160211          2539   CAA39089 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
7 ...

I am sure if you look at all the available filters and attributes, you can pinpoint the ID you are looking for.

Thanks. I can reproduce that, but only paint by numbers with R. How would I modify that to loop through a list of 4000 proteins and ensure the pdb ID was included in the output?

You can add 'pdb' to the list of attributes. That will add a column with the pdb id. c('pdb','hgnc_symbol',....),

Ok. Get it. Thanks.

I am not sure how your input file is (csv, vector, data.frame), but to run on multiple pdbs, you need to change yourvalues (values= c("2bhl","pdbID2","pdbID3",...))

My input is just a text file of filenames, one per line, generated from a unix ls command. I can convert it into a string of the type you suggest, although I wonder if there is any limit on the length of such strings. Referencing a text file would be preferable if it is possible.

I think the limit is 500 values x query, bu the newest versions of biomaRt should automatically split for you. I am not 100% sure, so you may want to double check. You can ls *.pdb >> all_pdbs.txt , use all_pdb_list <- read.delim("./all_pdbs.txt") to read all_pdbs.txt , and gsub to remove the extension from your entries (gsub(".pdb","", all_pdb_list)) .


text:
Melt p-values for CpG sites mapping to the same gene
I have some data I am working with, and I am curious if I am able to combine p-values from a paired t-test for CpG sites in the genome using Fisher's Method to get one p-value for each unique gene. Linked here is the Wikipedia page for Fisher's Method. I understand that an assumption of the method used is that each individual p-value being combined must be independent. I'm relatively new to biostatistics, so I'm curious if using CpGs from the same gene would violate this assumption.

Why do you want to melt CpG values? Do you want to know if an specific gene/region is more phosphorylated than others? What are your hypothesis (you tell about a t-test for the CpG sites? Are you testing if they are phosphorylated vs not or more than the rest of the genome or...?

In a nutshell, your intuition is correct: they are not independent.

Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.

What is typically done in methylation analysis is to assess the islands of methylations.
Check this workflow, in the section linked it uses some predefined islands for instance. I am no expert on this area but you could asses if the islands or certain regions are more methylated than expected.


text:
What is the best distribution to model the FPKM values from normalized RNA-Seq data?
I know that the discrete raw counts from the RNA-Seq data are usually modeled by a negative binomial or a Poisson distribution, but what I am working on are the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values which represent a transformation of the raw values. They are continuous except many actual 0s. What is the probability distribution FPKM values are likely to follow, or are assumed to follow? Even after taking logarithm twice (after replacing 0s with 1, of course), the data does not follow a normal distribution since there are a lot of 0s.

I think you may want to look into zero-inflated distributions?

Could you show or provide some sample data to explain with? What have you read about limma, DESeq2, DESeq ?

FPKM/TPM values are generally log-normally distributed. Reference : Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels

Sorry, I added a reference. I'm not sure why mRNA levels follow that distribution but it has been established. As they mention in the paper, it's a fairly common distribution in nature. On the other hand, 3' capture methods of single-cell RNA seq such as DropSeq usually yield a negative binomial distribution, presumably because of the high drop out rates of this method of transcript capture.

Well done. The answer looks WAY better with the reference. If you would like to push it even further. This type of information you have written in the comment you can add to answer as well to give a complete picture.

If some data follows a log-normal distribution, does it mean that it would follow normal distribution after you take the log of the values? That's not what I am seeing...

It should...what does expression look like for a gene you know is expressed in your sample after you log normalize your data?

Cuffdiff is one the tools that deals FPKM data for analysis of differential expression.By default they use the negative binomial with some condition specific parameter to avoid biases. There are some other models also available in Cuffdiff. Check the following documents for more details,
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html#library-types
https://www.reddit.com/r/bioinformatics/comments/3bx3em/fpkm_vs_raw_read_count_for_differential/

FPKM data is not composed of integers, it is interesting that they are using negative binomial distribution.


text:
SNP-phenotype association analysis only using the SNPs on a specific gene
I know that the GWAS association p-value threshold is 1e-8. This makes sense because in GWAS, you make ~1 million hypothesis (i.e. use that many SNPs in the association test). However, let's say I do not do genome-wide analysis but only use the 100 SNPs from a specific gene to check for association with my phenotype. Then, do you think a raw p-value of 1e-4 (which becomes 1e-2 after multiple hypothesis correction) is a reliable p-value to come to a conclusion that there is an association between that SNP and the checked phenotype? Obviously, if I checked all 1 million SNPs instead, this p-value would be insignificant after correction, but I am not checking 1 million SNPs, as I mentioned.

I know that the GWAS association p-value threshold is 1e-8

This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is unexpected shared genetic structure in the cases or controls, and this cannot be excluded by a p-value threshold, regardless of how low it is set.
In your 100-SNP example, you are correct that the bonferroni correction would be a threshold of 1e-4 (for a desired significance threshold of 0.01), but it's important to be cautious about GWAS results, even when the p-value indicates otherwise. Trust no one, especially yourself. P-values should never be used as an indication of importance (which suggests that most GWAS results should probably be reconsidered), and are better used in conjunction with other evidence for identifying relevant / significant tests.
My recommendation for GWAS (or any similar bulk genotyping study, as in your 100-SNP subset) is that at least a bootstrap sub-sampling process is carried out, to make sure that observations found within the comparison of interest are at least well replicated when comparing smaller sub-groups of the cases and controls. More information on that can be found in the poster I presented at Queenstown Research Week last year, and detail about the implementation in my draft preprint (which will probably never be properly published, given how old the research is).

Thanks for the reply, those make a lot of sense. A followup question: I was also told by a GWAS expert that when the minor allele frequency is low, one would not trust the p-value even if it is significant. He said so when I showed him the p-value of 1e-4 for a SNP with a 2% minor allele frequency. Do you think it is due to the need for the additional bootstrap analysis you suggested? I guess a low p-value would not be replicated in bootstrap analysis in case of rare variants and small overall sample size. Could you share your suggestions specifically for analyzing rare variants?

Whether or not it will be internally replicable will depend on the sample size of the subsampled population. It would be possible to specifically subsample to increase the frequency of single SNPs, but that might cause problems with associated interaction effects.

Certainly, you need to take more care about rare variants. Most tests will provide both a p-value and an effect size (which is better than the p-value for ranking SNPs). If the effect size is large, then it lends more weight to an association near the SNP locus. It may also be a good idea to do a haplotype association test within the haplotype block(s) that include the SNP, as any linked SNP could be causative. Bear in mind that GWAS are only testing for association, and there are other follow-up studies that should be done to support findings of association.


text:
Downloading CopyNumber450k package from Bioconductor
I am trying to download CopyNumber450k package in the R processor on my laptop in order to see if I could do a step of protocol for my research proffessor.  I was unable to successfully download the package I think (because it was too outdated for the newest version of R?) The version of R I have downloaded is 3.4.1. This package appears to work on the version 1.8.0 of R. I thought this was strange, as that seems really old. I then tried to see if this download was possible using some sort of e-version of an R-processor : https://rdrr.io/bioc/CopyNumber450k/   Even the e-version appeared to be unable to find and download this package of R. I have a feeling that I may be misinterpreting something? I tried to find an equivalent package that could run on the newer version of R and didn’t seem to find anything. Then, I tried to see if it was possible to download an old version of R. I couldn’t find a clear way to do that either. Does anyone on bioinformatics stack exchange have insight regarging what I may be doing wrong, something i may be misinterpreting, or how I may overcome this challenge?
Thanks

This is more of a troubleshooting question than a bioinformatics question. It would be better to search for your specific error on Google, than post a question here. You can probably find something on StackOverflow, or an R-specific site

I think you need to explain better what you have tried, including some code. You say that you were unsuccessful to download it, because of an outdated version of R. How did you do that?
I am not sure what you have tried, but if you read the instructions for installation here, you need to type the following in R console:
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("CopyNumber450k")

Show us what kind of error you get from this. Then we can maybe help you further.


text:
WGCNA co-expression network analysis with less than 20 samples
I have two cancer subtypes data. Subtype A is 14 samples and Subtype B is 23 samples. I'm interested in identifying the functions of some LncRNAs in the Subtype A group. For this I'm using all protein-coding genes and interested LncRNAs expression data (normalised counts).
From the WGCNA tutorial, I see that minimum number of samples needs to be 20. So, I guess using only Subtype A 14 samples for the analysis is not possible.
Can I use both Subtype A 14 samples and Subtype B 23 samples together for the analysis? In case if I use both, once I get the modules how do I know that the genes in modules are related to Subtype A or Subtype B? The pathways I get from the interested module is related to Subtype a or Subtype B? How this can be identified?

It's hard to say without knowing how different the subtypes are. If you do a common WGCNA, you may find modules related to the differences between A and B as well as to drivers of expression variation in A and B (common and separate). Since A has fewer samples, the variation will have to be stronger to be reflected by WGCNA modules.
You can run the WGCNA, find modules that vary the most in A samples (these may or may not also vary in B samples), then see whether LncRNAs belong to these modules and have a high kME (module membership).

thanks a lot Peter


text:
WGCNA: Problem with selecting soft threshold
I have 25 tumor samples with counts data. Initially, I filtered out low expressed genes and then converted counts to varianceStabilizingTransformation using deseq2 package. With this data I started using WGCNA for co-expression network analysis.
For selecting the soft threshold I see very strange plot. R2 cutoff is 0.8 and I see that none of the scale free topology model fit is above that.
Here is the code I used:
df is a dataframe with genes as rows and 25 samples as columns with counts data.
library("DESeq2")
filtered.counts <- df[rowSums(df==0)<5, ]

U3 <- as.matrix(filtered.counts)
vsd <- vst(U3, blind=FALSE)

oed <- vsd

gene.names=rownames(oed)
trans.oed=t(oed)
dim(trans.oed)

n=16462;
datExpr=trans.oed[,1:n]
dim(datExpr)

SubGeneNames=gene.names[1:n]

library(WGCNA)
options(stringsAsFactors = FALSE);
allowWGCNAThreads()

powers = c(c(1:10), seq(from = 12, to=20, by=2));
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
                      powerVector = powers,corFnc = cor,
                      corOptions = list(use = 'p'),networkType = "unsigned")

   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1 0.000273  0.041          0.786 4140.00  4000.000 6730.0
2      2 0.428000 -1.130          0.852 1570.00  1400.000 3770.0
3      3 0.673000 -1.540          0.882  729.00   579.000 2410.0
4      4 0.737000 -1.720          0.891  383.00   268.000 1670.0
5      5 0.745000 -1.830          0.886  220.00   134.000 1210.0
6      6 0.704000 -1.990          0.860  134.00    71.700  909.0
7      7 0.737000 -1.980          0.890   85.90    40.300  701.0
8      8 0.742000 -2.020          0.903   57.30    23.500  551.0
9      9 0.733000 -2.090          0.917   39.40    14.200  441.0
10    10 0.750000 -2.080          0.934   27.80     8.810  357.0
11    12 0.770000 -2.080          0.952   14.80     3.670  243.0
12    14 0.775000 -2.090          0.951    8.43     1.660  171.0
13    16 0.397000 -2.820          0.459    5.06     0.801  124.0
14    18 0.409000 -2.770          0.474    3.18     0.406   91.4
15    20 0.421000 -2.720          0.490    2.06     0.215   68.8

# Plot the results
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n", main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");

Is this the right way? Which soft threshold power should I select?
Any help is appreciated. thanq.

Don't worry about it too much. I would go with power 8 based on my general experience (also reflected in WGCNA FAQ) and on the mean connectivity around 50 and median around 20, which seems reasonable to me.

thanq very much Peter.


text:
How to deal with duplicate genes having different expression values?
I have RNA-Seq data which is FPKM. In the dataframe df first column is gene_name and the other 100 columns are samples.
Usually if it is counts data I do like following:
df2 <- aggregate(. ~ gene_name, data = df, max)

I'm not sure what do with the FPKM data if there are duplicate genes with different FPKM value for the same sample.
Lets say it looks like this:
gene_name     sample1        sample2        sample3
5S_rRNA      0.3206147    0.3327312      0.377578
5S_rRNA      0.3342000    0.0000000      0.1305166

Any suggestions please.

How were the read mapped? Do you have access to the raw counts data?

I have downloaded TCGA FPKM data. Would like to use the data for co-expression analysis with cemitool

I assume you're familiar with the various issues surrounding FPKMs, so I'll not expound upon them.
As a general rule, you should be using gene IDs rather than gene names, since the former are unique while the latter are not. If you only have access to data quantified on gene names, then the appropriate way to merge RPKMs is with a weighted sum:
$FPKM_{gene} = \frac{FPKM_{copy1} * Length_{copy1} + FPKM_{copy2} * Length_{copy2}}{Length_{copy1} + Length_{copy2}}$
As an aside, rRNA expression levels will tend to be wrong, since they're normally excluded by either poly-A enrichment or the use of ribo-zero.

Sure I will use gene ids. but lets say I have counts data and before doing differential analysis do you think this is the right to way to filter duplicates; aggregate(. ~ gene_name, data = df, max)

If you have counts then sum them.

So, if the first column is Ensembl ids and all other columns are samples in a dataframe "df", Is this right way to sum the duplicate Ensembl ids? aggregate(. ~Ensembl_ids, data=df, FUN=sum)

There should never be duplicate Ensembl IDs, they're unique. You only have issues like this with gene names or UCSC IDs.

oh yes, so with gene names Is this the way to sum the counts of duplicates aggregate(. ~gene_names, data=df, FUN=sum)

I have to say I disagree with the weighted approach - why would you place more emphasis on any of the individdual loci (in this case the longer)?

The best way to deal with this is to use unique gene IDs, for example ensembl accession numbers. So use the ensemble gtf annotation when quantifying the read counts and not the gene symbols. Just to illustrate, when I look for "5S_rRNA" in ensembl's annotation, i see 18 different "genes" with that gene symbol. But which 2 you have is unclear now.
grep "5S_rRNA" ensembl_symbol.txt
"ENSG00000252830"       "5S_rRNA"
"ENSG00000276442"       "5S_rRNA"
"ENSG00000274408"       "5S_rRNA"
"ENSG00000274059"       "5S_rRNA"
"ENSG00000276861"       "5S_rRNA"
"ENSG00000274759"       "5S_rRNA"
"ENSG00000280646"       "5S_rRNA"
"ENSG00000277411"       "5S_rRNA"
"ENSG00000201285"       "5S_rRNA"
"ENSG00000212595"       "5S_rRNA"
"ENSG00000277418"       "5S_rRNA"
"ENSG00000277049"       "5S_rRNA"
"ENSG00000274097"       "5S_rRNA"
"ENSG00000277488"       "5S_rRNA"
"ENSG00000274663"       "5S_rRNA"
"ENSG00000283433"       "5S_rRNA"
"ENSG00000275305"       "5S_rRNA"
"ENSG00000278457"       "5S_rRNA"

Yes you are right. But can I take median of those duplicates? something like aggregate(. ~ gene_name, data = df, median)

Expanding this answer, it seems that it was mapped to ensembl ids. So I would recommend to keep them as they are, and only when looking for biological meaning to change to gene symbols (or when talking to others)

I would strongly advice not to use FPKM at all, it is pretty clear now that they are not the right values to use for any analysis. If you still want to use it, keep your two genes apart like they were.


text:
How to get rid of the temp files while using "featureCounts" for extracting readcounts from bam files?
I used this command for extracting read counts for almost 60 samples using featureCounts.
featureCounts -a /gencode.v27.primary_assembly.annotation_nochr.gtf -F GTF -p -s 2 -T 8 -o counts.txt sample.sorted.bam

I see many "temp-core-100373-A4BF0116E148.sam-TH00-BK004414.tmp" files like this dumped in my directory. I'm able to delete them only manually. Is there a way to get rid of these tmp files from my directory while using the command?

You shouldn't delete them while featureCounts is running, it's creating read-name-sorted temporary files so it can more easily look at both mates together. Just wait for it to finish running and it will clean up the files itself. Doing so beforehand is liable to cause an error.

Thanks. But will there be any memory problem for those temp files? Asking this cz right now I have only 60 samples, for those samples I removed temp files manually. It didn't clean up even after finishing the job. I have to use featureCounts on at least 900 samples in future.

There aren't memory issues involved, though I suppose you could run into temporary space constraints. If featureCounts didn't clean up after itself then I encourage you to check and ensure it didn't encounter an error. I have never had it leave temp files around otherwise.


text:
Normalizing microarray data for clustering heat map
I wanted to generate a clustering heat map for the microarray data. This is the first time I'm working on Microarray data. I read some tutorials but have few doubts.
I'm using microarray (Affymetrix SNP 6.0 data) gene expression data. For example the data looks like following:
ProbeID S1  S2  S3  S4  S5  S6  S7
10008   131.4311    369.4926    222.0441    687.4181    176.8892    258.1233    316.5573
10010   78.73022    83.97501    81.56039    86.11443    78.09758    81.88231    84.17101
10014   90.02816    95.07267    101.1761    93.35585    81.96468    94.3553 93.89527
10017   79.86837    81.63064    88.19524    79.47265    76.3437 101.6351    93.71674
10019   99.03493    109.1835    104.97  102.7423    108.3677    98.93459    101.4052
10020   79.58075    84.28915    90.53562    74.47786    75.96112    96.39649    95.8828
10021   121.5373    149.9351    146.5956    122.8523    110.5759    132.4268    130.4409
10025   616.5994    1326.735    1358.187    2315.851    1068.745    3229.759    4435.021
10035   70.44073    69.56772    68.25446    68.35857    70.86771    74.3843 67.93569

I created expression matrix and did log2 conversion. Now the data looks like following.
            S1  S2  S3  S4  S5  S6  S7
Gene1   6.429276339 6.338451158 6.333760753 6.419191996 6.503471181 6.329103499 6.211373601
Gene2   6.379471993 6.296018518 6.237316465 6.2696332   6.329489132 6.359770303 6.240070336
Gene3   12.84498365 12.00265682 13.92741965 12.553162   13.39001307 13.8933423  12.58695704
Gene4   8.661860382 6.723004202 6.300975176 7.661012019 7.905219709 6.957578023 6.70945883
Gene5   6.945382967 6.814979733 6.575916303 6.63460198  6.627380524 6.733926424 6.280618235
Gene6   9.280222581 8.81560969  9.683073561 9.480038673 8.801438707 9.190943942 7.743705471
Gene7   6.593564871 6.63502488  6.389962535 6.511360029 6.694572404 6.6492763   6.527199544
Gene8   6.431615372 6.309515078 6.248153876 6.329288965 6.46768078  6.355268547 6.384284754
Gene9   7.513349406 7.234654595 7.490935892 6.801368215 7.323811386 7.018733196 7.055932044

Before clustering for normalisation:
I applied a function for normalisation
data = t(apply(data, 1, function(x) {
  q10 = quantile(x, 0.1)
  q90 = quantile(x, 0.9)
  x[x < q10] = q10
  x[x > q90] = q90
  scale(x)
}))

When I apply the above function I see there are values which are positive and also negative.
And when I used "rma" for background correction and "normalise" function from "oligo" package I don't see any negative values.
I don't understand.
Which one should I use to check the clustering among samples? Any help is appreciated.

You see negative values with your function because you're setting the average of each row to 0 and its standard deviation to 1.
In general, I would trust a standard normalization method (rma in this case) more than some random "truncate and then scale the rows" method. Your method isn't even doing any between-array normalization, which is the benefit of rma.

Thankyou for the reply. So, I will use the rma normalisation method.                                                 ty <- backgroundCorrect(data); normdata <- normalize(ty); And then using the normdata I can generate a clustering heat map. Am I right?

Could you please tell me whether what I mentioned in my previous comment is right or not?

That or rma() directly.

Ok. I did the normalisation like I mentioned in my previous comment. I also generated a clustering heat map. Could please have a look and tell me whether this is fine or not for clustering? https://imgur.com/3vZsO4Y

It looks like a clustered heatmap, whether it's "fine" is something that only you can answer.

Yes, I know what I mean is whether there is anything needs to be added for clustering heat map. Thank you !!

There isn't anything else absolutely required. It might be nice to not split by Type, since I imagine someone would like to see if the samples cluster according to it.

Ok. But when I don't use split it looks like this. https://imgur.com/a/MJotE

That's telling you something that's rather important to know.

This means samples are not clustered. Am I right?

They don't cluster by Type, at least.

Negative values are to be expected as you have log-transformed the data. Any data value below 2 will have a negative log2.
There is a difference between normalising for analysis and truncating for plotting. The RMA normalisation or similar is intended for differential expression analysis. I would recommend using this normalisation first.
When plotting, the normalised data may still contain outliers which skew the plotting parameters (such as the colour scale). Due to the large number of genes, there may be some far beyond 3 standard deviations and plotting these at their true value will make it hard to discern trends in the middle. Truncating the data for a heatmap is common, typically at |z| > 3 but quantiles should work much the same. This is a matter of visualisation, not clustering analysis. It would not be needed unless you were plotting the data.

Hello Tom, Yes I want to plot these data. You mean I should use the normalisation which is mentioned above in my question?

You should normalise the data first as you wish to analyse it further. The visualisations should as closely matched the data you've used, to adjust for batch effects, etc. If necessary to plot the values on a sensible range further standardisation (such as the z-transform) may be done. While there are best practices, expression analysis is not an exact science and involves a lot of decisions, depending on the variation you want to focus on, for what you can show from the data.