text:
What is deep sequencing?
People talk about deep sequencing. Is there any way to calculate how deep the sequencing is ? What should be the optimum depth to get reliable data ?
I am doing whole genome sequencing of a virus genome which is 10 kb long. I got 80000 reads from Illumina sequencing. Is there any way to tell about my sequencing depth ?

I found a post useful for this topic.
It explains the difference of coverage and depth.  It also has a useful explanation on how to calculate coverage and depth.

Here is a copy of what the link says just encase the post is removed:

Depth of coverage
How strong is a genome "covered" by sequenced fragments (short reads)?
Per-base coverage is the average number of times a base of a genome is
  sequenced. The coverage depth of a genome is calculated as the number
  of bases of all short reads that match a genome divided by the length
  of this genome. It is often expressed as 1X, 2X, 3X,... (1, 2, or, 3
  times coverage).
Breadth of coverage
How much of a genome is "covered" by short reads? Are there regions that are not covered, even not by a single read?
Breadth of coverage is the percentage of bases of a reference genome
  that are covered with a certain depth. For example: 90% of a genome is
  covered at 1X depth; and still 70% is covered at 5X depth.

SAMtools: get breadth of coverage
Get breadth of coverage of a reference genome.
How many bases of my reference genome are covered by my short reads sample at a level of 5X or higher?
 # set shell variables to your file-names

 REF_GENOME_FILE=my_ref_genome.fna
 WGS_SAMPLE=my_short_read_sample.fastq
 MIN_COVERAGE_DEPTH=5

Map short reads against reference genome (bowtie2)
# create bowtie2 index database (database name: refgenome)

bowtie2-build ${REF_GENOME_FILE} refgenome
ls
refgenome.1.bt2
refgenome.2.bt2
refgenome.3.bt2
refgenome.4.bt2
refgenome.rev.1.bt2
refgenome.rev.2.bt2

# map reads and sort bam file
bowtie2 -x refgenome --no-unal -U ${WGS_SAMPLE} -S - -p 12 | \
samtools view -bS - | \
samtools sort -m 5G - mapping_result_sorted.bam

-U ${WGS_SAMPLE} consider all reads as unpaired reads
-S -             bowtie output in SAM format, write to standard out
-x refgenome     use reference genome bowtie2-database
--no-unal        suppress SAM records for reads that failed to align
-p 12            use 12 processors

Get breadth of coverage (SAMtools)
# create samtools index
samtools index mapping_result_sorted.bam

# get total number of bases covered at MIN_COVERAGE_DEPTH or higher
samtools mpileup mapping_result_sorted.bam | awk -v
X="${MIN_COVERAGE_DEPTH}" '$4>=X' | wc -l
32876

# get length of reference genome
bowtie2-inspect -s refgenome | awk '{ FS = "\t" } ; BEGIN{L=0}; {L=L+$3};
END{print L}'
45678

Result: breadth of reference genome coverage
total number of covered bases: 32876  (with >= 5X coverage depth)
Depth of coverage (average per-base coverage): 0.719 X (32876 ÷ 45678) (total number of covered bases divided by reference genome length)
percent: 71.9%  (0.719 × 100) (coverage breadth in percent)

Nice picture! Could you include the relevant parts of the linked page in the answer to make the answer self sustained (and avoid having to go one link further).

Thanks for the answer Samantha. This exactly answers my question.

There are several questions in your post I'll try to answer each one:

Is there any way to calculate how deep the sequencing is ?

See gringer's answer. TLDR: The depth of the sequencing is how many times each position has been sequenced.

What should be the optimum depth to get reliable data ?

The optimum depth depends on what you want to do with that data. Usually,  30X or 50X is done, but in targeted sequencing you can sequence a lot more reaching 1000X (which would be really deep).

Is there any way to tell about my sequencing depth ?

If you know the length of the region sequenced and how many bases you have in the reads, you can calculate the general sequencing depth.

Sequencing depth is typically calculated as the number of total bases sequenced divided by the number of bases in the target genome. An Illumina sequencing run with 2x125 bp reads and 500 million read pairs sequenced would be a sequencing depth of about 40X (assuming my calculations are correct for a 3 billion base genome).
The sequencing depth depends on how long each read is, as well as whether it is a paired-end sequencing run. 80,000 read pairs with each having 125 bp would be a coverage of 2000X for a 10 kbp viral genome (assuming the sample were entirely composed of that virus).
The optimum depth changes depending on the application. A genotyping-by-sequencing analysis, combined with imputation, could work well with 96 samples at a 1X coverage (or possibly less). On the other end of the scale, a cell-free DNA analysis to identify structural variation in metastasing cancer may need a few thousand times coverage of the host DNA.


text:
Adding entries to bigwig file
I generate bigwig files using a shell script based on bedtools genomecov (to generate a bedgraph from a bam file) bedmap (to compute means across 10 bp bins) and bedGraphToBigWig to convert the binned bedgraph into bigwig.
Sometimes, the bam file has no data on some chromosomes, this results in a bedgraph file with no entries for the empty chromosomes.
I would like to add the missing entries. How could I do this?
I tried to see what I could do with pyBigWig, opening the bedgraph file in "a" mode, but don't really know how to proceed:
# I have the chromosome sizes for my genome of interest already in a list:
chrom_sizes
#[('I', 15072434),
# ('II', 15279421),
# ('III', 13783801),
# ('IV', 17493829),
# ('MtDNA', 13794),
# ('V', 20924180),
# ('X', 17718942)]
bw_test = pyBigWig.open("path/to/bigwig", "a")
bw_test.chroms()
# {'IV': 17493829, 'V': 20924180, 'X': 17718942}
bw_test.addHeader(chrom_sizes)
# RuntimeError: Received an error in bwCreateHdr

Presuming you still have your BAM file sitting around somewhere:
#!/usr/bin/env python
import pysam
import pyBigWig

bam = pysam.AlignmentFile("alignments.bam")
ifile = pyBigWig.open("your current file.bw")
ofile = pyBigWig.open("your new file.bw", "w")

# Create/Add the header
hdr = []
for chrom in bam.header['SQ']:
    hdr.append((chrom['SN'], chrom['LN']))
bam.close()
ofile.addHeader(hdr)

# Write the bigWig
for _ in hdr:
    try:
        ints = ifile.intervals(_[0])
        starts = [x[0] for x in ints]
        ends = [x[1] for x in ints]
        vals = [x[2] for x in ints]
        ofile.addEntries([_[0]] * len(starts), starts, ends=ends, values=vals)
    except:
        # intervals() will throw an error if the chromosome isn't present
        pass
ofile.close()
ifile.close()

The reason you can't append to a bigWig file is that it has a various indices and zoom levels typically built in. You'd have to add to them without overwriting neighboring blocks of data. In theory this might sometimes be theoretically possible, but not typically and I'm not aware of tools that support this.

If I understand correctly, the bam file is not necessary in my case, because I already have the chrom_sizes list, which should be the same as the hdr you obtain by parsing the bam header.

Correct, if you already have such a list then you don’t need a BAM file.

Just a detail: It seems strange to me to use a _ placeholder as looping variable if you then use it to access its first element. I suppose for (chrom, _) in hdr: would work as well while avoiding you the _[0].

Yup that’d work too.

source: Stack Exchange

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

sort: votes


text:
Differential Expression With 2 Treatments
I have a expression data from a small cohort of samples taken at baseline and after 2 independent treatments. I can do differential expression contrasting T1 and T2 or I can contrast T1 vs baseline and T2 vs baseline and look at the differences. What is the difference between these two analyses? Are they both valid? If so what inferences can be drawn from each?

If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline.
The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes differentially expressed vs. baseline in at least one condition (e.g., to get rid of genes only slightly higher due to T1 and slightly lower in T2, but not different enough in either case to be DE). Use a pretty lax p-value threshold (e.g., 0.1 or 0.2) for this.

Thanks. What is your opinion on genes that may be significant in just 1 of the comparisons to baseline but not significant in the direct T1 vs T2 contrast?

It there's no difference in the T1 vs. T2 contrast then there's no difference. It doesn't matter if you happen to get a significant difference in one of the conditions vs. baseline. You only want to use the comparisons to baseline for filtering your T1 vs. T2 results (assuming you have enough DE genes to bother doing so).

One more question. So if I wanted to do a pathway analysis and compare the two treatments that way. When I generate my gene lists how should I treat genes that are DE in only one treatment vs baseline but not DE in the direct T1 vs T2 contrast?

Treat them as not DE.

So any genes that have similar responses from baseline in both treatments get discarded? Wouldn't you be missing a lot of context if a gene has to be DE vs baseline in both treatments and also DE directly with each other to make either final list?

Yes, any genes with similar responses should get discarded, since they're not different. A gene doesn't have to be DE in either treatment vs. baseline to be DE between them. That's also why I suggested using a very lax p-value threshold if you want to do any filtering of the T2 vs. T1 results.

Sorry, I meant doing a pathway analysis of each treatment vs baseline separately to find the raw response of each treatment. Since they are different from baseline I would still include them in only their respective pathway analysis?

Yes, certainly include them in such cases.

Cool, thanks a lot for your help. I really appreciate it.

Depending upon the type of treatment used the set of DEGs will change.If the treatments have similar kind of effect you will get a small list (less variable genes will have higher p-val) using a cutoff of p<0.05. So, it's better to start with control vs treated comparisons then T1 vs T2.
Compare the lists form CTRL vs T1 and CTRL vs T2 you will get the genes that are expressed in both conditions as well as unique to individual treatments.

If a gene has a p-value of 0.049 in CTRL vs. T1 and 0.05 in CTRL vs. T2 you would be categorizing it as uniquely DE in one condition, when there's no actual difference between the conditions.

Yes, it's a tricky situation.

I agree it is a trick situation in my case it is modes of exercise so they are extremely similar responses. The T1 vs T2 comparison yields almost nothing where the vs baseline comparisons are quite different.

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


text:
How to prevent sklearn Imputer(missing_values="NaN", strategy="mean", axis=0).fit_transform(data) from removing columns with only NA in them
I am trying to test a preexisting python machine learning script with a subset of my genetic data.  One of the feature columns I am using happens to only have NA values in it.  I lose this column when I use the following code:
Imputer(missing_values='NaN', strategy='mean', axis=0).fit_transform(data)
This causes problems later on because my model expects the correct number of features, and I do not want to retrain my model at this time.  Is there any way to work with Imputer or fit_transform to prevent this column from getting dropped?

Join the "NA" containing column after the transformation.

I am in no way expert in machine-learning but this seems more like a programming question than a bioinformatician problem, (although you use genetic data). Consider "telling" the model that there is now one feature less. But if you have new data enough different from the original you need to retrain your model AFAIK...

We might conceivably be able to help if you gave us some details but without your code or a dataset to test it on, we can't. You are using various commands there which, presumably, are defined in your model. We can't help if you don't tell us how they work or what they are.

authors: أماره أبو صهيون


text:
Coverage calculation: long reads (RNA-seq)
Say your aim is to calculate the coverage of an RNA-seq experiment generated with long-read sequencing (so, uneven read length).
Up to now, I relied on the Lander/Waterman equation:
$$C = L*N / G$$
where:
$C$ = final coverage
$G$ = haploid genome (or transcriptome) length
$L$ = read length
$N$ = number of reads
I have two major conceptual issues with this formula:

$L$ is uneven in this case. Should we use the median/mean read length instead? If so, side question: if the RNA-seq experiment was executed following the size-fractionated protocol (PacBio), is it advisable to merge all the fractions prior to coverage calculation or should it be calculated separately for each fraction?
Should $N$ consist of the mapped reads to the reference or simply the number of reads as found in the subreads whole dataset?

The read length is irrelevant when calculating the mean coverage statistic. It's simply the total number of bases sequenced divided by the target Xome length.
In the example provided in the question, $L*N$ is equivalently expressed as $\sum^{N} l_x$, or the total length of sequenced reads. For equal-length reads, this is easier to calculate as number of reads $*$ length of reads, but the full sum is necessary if read lengths are not equal.
For what it's worth, this value, $L*N$, is also equivalent to the average read length multiplied by the number of reads, i.e.
$$\frac{\sum^{N} l_x}{N} * N$$
... which makes it a bit more obvious why the total length of sequenced bases is very slightly easier to calculate.
In answer to the second part of the question, $N$ should ideally correspond to only the sequences that are from the reference. For an assembly, only the reads that were used to create the assembly should count when considering coverage. When there's no likelihood of contamination (and especially for RNASeq with variable coverage per transcript), a statistic using the total number of reads should be good enough for most purposes that I can think of. High-throughput sequencing is such a fuzzy process that a few tens of millions of bases here or there shouldn't make a difference.

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


text:
What is this 5-column sequencing format from 2009?
There is high throughput sequencing data here, and I don't know what format it is in.
It was submitted in 2009, and the description says the following:

Library strategy:     ncRNA-Seq
Library source:   transcriptomic
Library selection:    size fractionation
Instrument model:     Illumina Genome Analyzer II
Description:  CIPPNK, tar file of Illumina *_seq.txt files provided as supplementary file

I got the archive here:
ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM455nnn/GSM455387/suppl/GSM455387%5FWT%5FCIPPNK%5Fseq%5Fs1%2Etar%2Egz
Inside, there are 330 files from s_1_0001_seq.txt to s_1_0330_seq.txt which are tab-delimited text files where the first column is always 1, the second has the number found in the file name, then 2 mysterious integers, and then what looks like a read of length 36, with sometimes a dot instead of a letter:
$ head s_1_0330_seq.txt
1   330 690 785 TTCCTACATTGTTCCCCCATGCTGTTGGCACCATCA
1   330 44  145 TTTTTATCACGAGTTTTAAATCTGTAGTCACCATCA
1   330 53  141 AATAATGCATAACAAAACGGAATCTGTAGAA.AAA.
1   330 784 461 TAATTGTAGTGATTGATCAATCTGTAGGCACCATCA
1   330 588 634 TATTATGCACATTTTCTAGTTCACTGTAGGCACCAT
1   330 718 678 TTACATGTTTCGGGTAGGAGCCTGTAGGCACCATCA
1   330 635 834 TGTGATCATTAGTTCAAAGCCCCCTGTCGGCACCCT
1   330 494 523 TGAAAATCAAAAATGCTGAACTGTAGGCACCATCAA
1   330 393 783 TTTTTTTTTAAATTTAAAAAAACTGTAGGCACCATC
1   330 48  148 GTTTAACCGTGTAGACGTTGGTTTCTGTAGGCACCA

This seems to be some "old" high-throughput sequencing format. I think someone in a 2008 message in sequanswers was dealing with this type of file:
http://seqanswers.com/forums/showpost.php?p=1841&postcount=8
What is this format that seemed so standard back then, that the authors did not give more info than describing the files as "Illumina *_seq.txt files"? I do not dare ask them for such a trivial question (the given contact is a Nobel prize laureate and is probably too busy to answer random bioinformatics questions).
In particular, what are the columns 3 and 4, and what do the dots mean?

This is an early Solexa/Illumina sequencer format. The columns are the identifier to location on the flowcell.  I believe the "." was the original placeholder for an unread base, which has been replaced by and "N" in current Illumina sequencing output.
From on http://www.crg.eu/en/content/processing-and-analysis-illumina-sequencing-data.

seq.txt (Gerald)
4 1 23 1566 ACCGCTCTCGTGCTCGTCGCTGCGTTGAGGCTTGCG
The first 4 columns are an ID that provides location details of the cluster in the flow cell, followed by the read sequence. Seq.txt file size typically ranges between 250 Megabytes and 1 Gigabyte. One file is generated per lane.

Ah, so if the authors did not provide other files, I guess the only way I can convert this to fastq is to chose arbitrary qualities.

Correct this format did not contain quality scores.  There was a qseq.txt file that did.  Although to be honest back then the quality scores were no that meaningful as the Illumina algorithm was not very accurate in assigning them.


text:
PLINK clump behavior on missing SNPs?
I have a long list of autoimmune-associated SNPs, and I want to boil it down so that I get one SNP representing each LD block. I chose to use PLINK's --clump option for this. I'm roughly following this tutorial (but analyzing my own data).
I don't have the full data from the original GWAS's, so I am using 1000 Genomes data as the reference for LD estimation. But, the 1000 Genomes data don't have all the SNPs I need. I end up a few hundred warnings:
Warning: 'rs12345678' is missing from the main dataset, and is a top variant.
How does PLINK handle this situation? More importantly, how should I handle this situation? I don't want to omit those SNPs, but I also don't want to take them all, for the same reason that I am doing this in the first place.

Are you sure they're missing? I mean, are you sure it isn't just a question of wrong/outdated/changed rsID? I've never done anything like this, but isn't there a way of using actual SNPs (genomic coordinates, reference allele and alternate allele) somehow? Perhaps by feeding it a vcf file?

That's a good point. I didn't realize rsID's could go stale like that.

@terdon On the contrary. rsIDs are designed to be a stable alternative to just using position and allele, so that when reference assembly changes, SNPs can be tracked by rsIDs.

@juod yes but, alas, they can change. Certainly in cases where rsIDs are merged for example. In those cases, one of the merged IDs is rendered obsolete.

Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's small or from a different population. The estimate will be biased, but you don't really need high precision for clumping.
Then again, you can skip PLINK and try searching for SNP prioritization methods based on summary statistics - it's an area of active development, and something might be available to suit your needs.

Unfortunately, I do not have even an individual-level dataset. But, I am accepting this answer because so many SNPs are missing, and I doubt different dbSNP versions can account for so many missing SNPs.

You can use imputation for guessing at what the missing SNPs might be based on known LD patterns in populations. This procedure will give you an idea of whether the recombinational history of genomic regions is important, rather than specific SNPs. Given that you're looking at LD-based measures, this should be okay (although the logic behind the calculation is a bit circular).
If you don't have a big compute cluster available, the Michigan Imputation Server is a free imputation service that can be used for academic purposes. Data is fed into the service as bgzip-compressed VCF files (output from plink as an uncompressed VCF file, then use bgzip to compress the files), and returned as a larger dataset with imputed SNPs derived from the HRC or 1kg reference panels. There's a bit of initial legwork required: you need to make sure that the variants as defined in the VCF file have the same location and strand as defined in the reference panel, otherwise you'll wait weeks in the queue and end up empty handed. To avoid hassle, I'd recommend removing any complementary SNPs (e.g. A/T, C/G), then flipping the strand for any remaining genotypes where the alleles don't match up.
I've recently done this again. Here's a summary of the things required to get MIS imputation working properly from a PLINK dataset:

Retrieve site data for HRC panel
Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
extract out variant names (remove duplicated rs numbers)
change 'RS' to 'rs' in plink BIM file
filter PLINK file against variants
adjust position of variants to fit with HRC panel
sort HRC variants by rs number
identify variants that need flipping
flip variants using PLINK, create combined VCF file
create VCF index for bcftools
extract per-chromosome VCF files
convert chromosome '23' back to chromosome 'X'
upload to MIS, remembering to separate out autosomal and X-chromosomal files

And a full example:
## Retrieve site data for HRC panel
wget 'ftp://ngs.sanger.ac.uk/production/hrc/HRC.r1-1/HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz'
## Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
zgrep -v '^#' HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz | cut -f 1-5 | \
  awk '{if(!($4 ~ /[AT]/ && $5 ~ /[AT]/) && !($4 ~ /[CG]/ && $5 ~ /[CG]/)){print}}' | perl -pe 's/^X/23/' > HRC_complementary_variants.tsv
## extract out variant names (remove duplicated rs numbers)
cut -f 3 HRC_complementary_variants.tsv | grep '^rs' | uniq -u > HRC_complementary_variants_namesOnly.txt
## change 'RS' to 'rs' in plink BIM file
perl -i -pe 's/RS/rs/' plinkFileHg19.bim
## filter PLINK file against variants
plink --bfile plinkFileHg19 --extract HRC_complementary_variants_namesOnly.txt --make-bed --out noComp_plinkFileHg19
## adjust position of variants to fit with HRC panel
plink --bfile noComp_plinkFileHg19 --update-chr HRC_complementary_variants.tsv 1 3 \
  --update-map HRC_complementary_variants.tsv 2 3 --make-bed --out repositioned_nc_plinkFileHg19
## sort HRC variants by rs number
sort -k 3,3 HRC_complementary_variants.tsv > rsSorted_HRC_complementary_variants.tsv
## identify variants that need flipping
join -1 3 -2 2 <(pv rsSorted_HRC_complementary_variants.tsv) <(sort -k 2,2 repositioned_nc_plinkFileHg19.bim) | \
  awk '{if(!($4 == $9) && !($4 == $8)){print $1}}' > toFlip_HRC_complementary_variants.txt
## flip variants using PLINK, create combined VCF file
plink --bfile repositioned_nc_plinkFileHg19 --flip toFlip_HRC_complementary_variants.txt --recode vcf bgz --out flipped_rnc_plinkFileHg19
## create VCF index for bcftools
bcftools index flipped_rnc_plinkFileHg19.vcf.gz
## extract per-chromosome VCF files
for x in $(seq 1 23); do echo ${x};
  bcftools filter -r ${x} flipped_rnc_plinkFileHg19.vcf.gz -O z -o chr${x}_flipped_rnc_plinkFileHg19.vcf.gz;
done
## convert chromosome '23' back to chromosome 'X'
zcat chr23_flipped_rnc_plinkFileHg19.vcf.gz | perl -pe 's/ID=23/ID=X/;s/^23\t/X\t/' > chrX_flipped_rnc_plinkFileHg19.vcf
bgzip chrX_flipped_rnc_plinkFileHg19.vcf

fyi, plink can directly generate bgzipped VCF files (--recode vcf bgz), and it'll use multiple threads to do the compression.

Yes, but the plink-generated bgz-VCF files didn't work properly with the Michigan Imputation Server last time I tried, and it's a real pain testing that out due to the queue length.

Thanks for mentioning this, I'll investigate.  I noticed that the recommended pipelines now include a "bcftools sort" step after VCF file generation even though that VCF would be position-sorted, so there may be a poorly documented assumption about how multiple variants sharing the same bp coordinate must be sorted, and/or an assumption about chromosome order.


text:
Somatic tumor only variant calling?
I'm evaluating possibilites for somatic tumor variant calling without paired-normal samples. I'm aware of the consequences without a normal sample.
All the popular variant calling tools such as Strelka, VarScan etc require a normal sample.
What are the some of the possibilites I have? I'm looking for good sensitivity in exchange for higher false positives.

There is a recent paper that attempts to do this:

ISOWN: accurate somatic mutation identification in the absence of
  normal tissue controls.
In this work, we describe the development, implementation, and
  validation of ISOWN, an accurate algorithm for predicting somatic
  mutations in cancer tissues in the absence of matching normal tissues.
  ISOWN is available as Open Source under Apache License 2.0 from
  https://github.com/ikalatskaya/ISOWN

I've been using the LoFreq* caller for exactly this. It is designed to find variants with very low frequency, so is well suited for this type of analysis.

LoFreq* (i.e. LoFreq version 2) is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of errors inherent in sequencing (e.g. mapping or base/indel alignment uncertainty), which are usually ignored by other methods or only used for filtering.
LoFreq* is very sensitive; most notably, it is able to predict variants below the average base-call quality (i.e. sequencing error rate). Each variant call is assigned a p-value which allows for rigorous false positive control. Even though it uses no approximations or heuristics, it is very efficient due to several runtime optimizations and also provides a (pseudo-)parallel implementation.

(In my experience it is not fast at all, but that's what they claim).
I have used it many times to call somatic variants using the germline caller (and not its somatic one which requires a tumor/normal pair). Because the tool is so sensitive, its germline caller actually deals quite well with somatic samples.
You can find instructions for its use here.

I have used LoFreq as well and think it does a good job.  It is very sensitive and almost always requires some post filtering to get a reasonable number of calls.  One hazard I have found with its model is that if the variant is actually the majority call at the position it can have trouble reporting it.  This can happen with samples that have a large tumor cellularity.

@Bioathlete oh? Thanks, that's good to know!

All the popular variant calling tools such as Strelka, VarScan etc
  require a normal sample.

Strelka and VarScan require a normal sample in somatic mode, but they both have a germline mode for unpaired analysis.
One of the most popular somatic variant callers is MuTect and that does not require a normal.

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


text:
Retrieve abstract/summary for NCBI bookshelf via Entrez
I want to retrieve abstracts/summaries for NCBI bookshelf entries, eg: "NBK1440"
Docs say the dbname is "books" and efetch guide says a rettype of "docsum" works for all databases. However when I make the query:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&rettype=docsum

returns:
<ERROR>UID=1440: cannot get document summary</ERROR>

The default query is useless (just echos back the ID):
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&retmode=text

returns
1440

Picking another 3:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1391,NBK410087,NBK1391&rettype=docsum

Fails the same way:
<ERROR>UID=1391: cannot get document summary</ERROR>
<ERROR>UID=410087: cannot get document summary</ERROR>
<ERROR>UID=1391: cannot get document summary</ERROR>

The same goes for esummary:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=1440

Dies on NBK1440 and removing NBK gives:
<ERROR>UID=1440: cannot get document summary</ERROR>

The values come from the (9Mb) ClinVar citations file:  https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/var_citations.txt - the first is from tihs line:
15048   9   1800562     NCBIBookShelf   NBK1440

ANSWER
Thanks to arupgsh - yes the ClinVar citation_id was a RID for bookshelf not an ID like it is for Pubmed and PMC. FYI my solution was:
def get_summary_for_bookshelf_rid(bookshelf_rid):
    handle = Entrez.esearch(db="books", term=bookshelf_rid)
    search_results = Entrez.read(handle)
    handle = Entrez.esummary(db="books", id=','.join(search_results['IdList']))
    results = Entrez.read(handle)

    for r in results:
        if r["RID"] == bookshelf_rid:
            return r

    return None

get_summary_for_bookshelf_rid("NBK1440")

returns:
{'PubDate': '2000/04/03 00:00', 'Title': 'HFE-Associated Hereditary Hemochromatosis', 'Text': '', 'RType': 'chapter', u'Id': '1475938', u'Item': [], 'Book': 'gene', 'Parents': '', 'RID': 'NBK1440', 'ID': 'gene/chapter/hemochromatosis/PMC'}

But unlike PubMed article, there is no entry for abstracts for NCBI bookshelf.

Entrez requires the unique identifiers(UID) for fetching related info. The id you are using in the query is RID. I guess that why you are getting <ERROR>UID=1440: cannot get document summary</ERROR>.
List of example uids using the term "ncbi+blast":
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=books&term=ncbi+blast
Summary of the first item:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=4459608
Using the Entrez python package:
from Bio import Entrez
Entrez.email = "youremail@domain.com"
handle = Entrez.esummary(db="books", id="4459608,4456634", retmode="xml")
records = Entrez.parse(handle)
for record in records:
     print(record['Title'])
handle.close()

I'm not sure about the availability of book abstracts.

Thanks - yes the ClinVar citation_id was a RID for bookshelf not an ID like it is for Pubmed and PMC

As arupgsh says, you need to use esearch to get a list of unique identifiers before using efetch to retrieve info about each result. I think the easiest way to do this is to use Entrez Direct, which allows you to simply pipe esearch output to efetch:
esearch -db books -query NBK1440 | efetch -format docsum

or
esearch -db books -query NBK1440 | esummary


text:
Biopython code refactoring
I want to make function get_abstract, but i don't know what this can look like if the function must return idlist.
from Bio import Entrez, Medline, SeqIO

list_of_ids = []
Entrez.email = 'ski89@gmail.com'

handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retmax=100000)
record = Entrez.read(handle_one)
idlist = record["IdList"]
for ID in idlist:
    list_of_ids.append(ID)

    if len(list_of_ids) == 100000:
        handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=100000, retmax=100000)
        record = Entrez.read(handle_one)
        idlist = record["IdList"]
        for ID in idlist:
            list_of_ids.append(ID)

            if len(list_of_ids) == 200000:
                handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=200000, retmax=100000)
                record = Entrez.read(handle_one)
                idlist = record["IdList"]
                for ID in idlist:
                    list_of_ids.append(ID)

so the args will be the start and TERM:
get_abstracts(start,TERM):
handle_one = Entrez.esearch(db="pubmed", term=TERM,
                               retmode="text", retstart=start, retmax=100000)
                    record = Entrez.read(handle_one)
                    idlist = record["IdList"]
                    for ID in idlist:
                        list_of_ids.append(ID)

Can someone help refactor this? cheers

Do you really mean to have a for loop nested inside another for loop that both reuse the same variables?

This type of questions are more suited for code review SE. What is your bioinformatic question? Do you need a new method to do it? A faster one?

Refactoring is for speed up script. I dont know any alternatives for download  many abstracts. Selenium is slow. And wanna be better programmer , thats why I ask. I dont have any helper.

Here is a sample of code that should get you started.
I created a separate get_abstracts function.  Then that function get called repetitively until less than 100K ids are returned.  If you wanted to get fancy this whole process could be called as a recursive function but that seemed overly complex.
BUFFER_SIZE = 100000

def get_abstracts(start, term):
    """Get N number of records from Entrez
    :returns list: List of ids returned from search
    """
    handle_one = Entrez.esearch(db="pubmed", term=term,
                       retmode="text", retstart=start, retmax=BUFFER_SIZE)
    record = Entrez.read(handle_one)
    idlist = record["IdList"]
    return idlist

list_of_ids = []
start_at = 0

# Iterate the search returning 100000
while true:
    ids = get_abstracts(start_at, 'e coli')

    list_of_ids.extend(ids)

    # Break if less than 100K results were returned
    if len(ids) < BUFFER_SIZE:
        break
    start_at += BUFFER_SIZE

# Do something with list_of_ids


text:
How to release an R package with several bed files?
I'm currently creating an R package, and my scripts require before any analysis is done that 1-2 bed files be loaded. Normally, I would run the scripts with the following:
library(data.table)
session_bed = read.table("/path/name/to/bedfile.bed")

and then refer to session_bed throughout the analysis.
However, in order to release this as an R package, I believe that these files must be downloaded at the time on the package installation.
(1) Normally, how are bed files distributed in R packages? Are these gzipped? Where are they placed?
According to http://r-pkgs.had.co.nz/data.html, all .bed files should be inst/extdata I guess. However, it's not clear to me how to "link" to this data within the R package.
How have other packages dealt with this? (Is there a size limit?)
(2) In this case, should the function loading the bed file be executed about loading the R library, library(package_name)? Is there more appropriate option?

You can store it in data/ subfolder. If one of your functions need this bed file, you can include importing it in your function.
Here are some examples how to get the data from data/ subfolder.
There are probably more ways to do this, this is just one way.

As benn said you can store them in data/ subfolder of your package. There is no size limit (AFAIK) but if you plan to publish your package in some repositories they might have limits or other restrictions. Bioconductor requires to divide the data required for a software package at certain threshold, but for small amount of data it is fine.
You can find them via system.file like: read.csv(system.file("data/mydata.csv", package = "mypackage")).
Of course your bed files will need the apropiate function to load them.
You can load when your library is called by modifying the local variable, or using .First (I am not so sure of the last one), but this pollutes the environment and sometimes if preferably to use data function. You might need to convert your files to .Rdata file.


text:
Simple Affy Differential Expression
I think I have a fairly easy-to-solve problem.
I have an expression table generated from an Affy set. It has one row of gene names (already symbols) and then 180 sample rows. 181 rows total.
So I made another excel file with one row with three different phenotypes I want to compare. That row has a label ("classes") and 180 names. 181 rows total.
In R, I can import each of these
dat <- read.csv("expression.csv")
dat <- data.matrix(dat)
pheno <- read.csv("pheno.csv")
pheno <- data.matrix(pheno)

I can then fit them
fit <- lmFit(dat,pheno)

And ask for an eBayes from limma
fit <- eBayes(fit)

Now, when I view the topTable,
table <- topTable(fit)

My first column has numbers (like "834","1142"), when I'd expect to have gene names there (like "ACTR", "CHRNA1").
What am I doing wrong? Or, rather, am I doing anything right?
Thanks!

What are your rownames of pheno? BTW are you sure that pheno has only the phenodata relevant for your differential expression analysis? Which is your model?

What do the first few lines of your lmFit input variables look like (i.e. dat and pheno)? Are the rows numbered, or do they use the actual symbols?
By default, I don't think read.csv expects a file to have row numbers (even though write.csv puts them in by default). Try adding a row.names parameter to the expression matrix read call:
dat <- read.csv("expression.csv", row.names=1)


text:
sum of products of mapq and mapped bases for each read in a from a BAM file
Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string).
For example, given two reads like this:
read1: mapq=40, mapped_bases=10
read2: mapq=20, mapped_bases=50

The number I'm looking for is:
(40*10) + (20*50) = 1400

One option here is to stream the file and get the two quantities line-by-line. But this seems likely to be slow. Thoughts welcome.

BAM files are stored line by line, so your only option is to process it that way. In python with pysam:
import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
        continue
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
bam.close()
print("the total is {}".format(tot))

One possibility uses bamutils (http://ngsutils.org/modules/bamutils/removeclipping/) to remove the clipped bases, then awk to get the sum of the products, something like:
bamutils removeclipping val.bam val_rc.bam
samtools view val_rc.bam | awk '{print length($10)*$5}' | awk '{sum+=$0} END { print "mapq*len sum =",sum}'

the bamutils step is rather slow, but it will do the job.

Why use three separate awk calls instead of samtools view val_rc.bam | awk '{sum+=length($10)*$5} END { print "mapq*len sum =",sum}' ?


text:
How to resolve in snakemake error : "Target rules may not contain wildcards."
I would like to do easily reproducible analysis using publicly available data from NCBI, so I have chosen a snakemake.
I would like to write a single rule, that would be able to download any genome given a species code name and separated table of species and their NCBI IDs. So I wrote a script scripts/download.sh that takes a species code and download the genome to data/<species_code>/genome.fa.gz. The script internally reads the table tables/download_table.tsv, where are corresponding species code names and NCBI IDs.
So I tried to do a snakemake like this :
species='Cbir Avag Fcan Lcla Dcor Dpac Pdav Psp62 Psp79 Minc1 Minc2 Mjav Mare Mflo Mhap Pant'

rule download:
    input :
        "tables/download_table.tsv"
    output :
        "data/{sp}/genome.fa.gz"
    shell :
        "scripts/download.sh {sp}"

However, snakemake returned an error message I do not really understand :
Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards.`.

Is there a way to write a single rule for downloading all the genomes?

The problem is that you need a master rule that requires all of your desired outputs as inputs, in your case it would be :
rule all:
    input:
        expand("data/{sp}/genome.fa.gz", sp=species.split(' '))

You'll also need separate download link inputs for each species. You could make a separate download_table.tsv for each species, but it would probably be easier to make a config file with this information, and add a params keyword to your rule. Something like:
rule download:
   params:
       url=config['locations']['sp']
   output :
       "data/{sp}/genome.fa.gz"
   shell :
       "scripts/download.sh {params.url}"

I wouldn't erase it if I were you. A lot of people get confused about Snakemake error messages, so its helpful to have some examples on here.

Then I will at least dramatically edit it. (I try to do it in a way so your answer still fits)

@KamilSJaron I agree the question should be kept, but maybe change the title to reflect the error message, this has nothing to do with NCBI

Also note that Snakemake has an NCBI remote provider: http://snakemake.readthedocs.io/en/stable/snakefiles/remote_files.html#genbank-ncbi-entrez

@JohannesKöster I was looking at it. But since I could not make it work without remote I have not even tried with remote. And then the problem I originally had with non-remote solution is very much unrelated to downloading of sequences. It's a good tip though, maybe work an answer?

Given the Makefile you provided and which someone then deleted from the question, you should probably just add this line (this should be the first rule):
rule all: input: ["data/{}/genome.fa.gz".format(x) for x in species.split()]

This rule just specifies a list of expected output files and corresponds to the following lines from the original Makefile:
GENOMES=$(patsubst %, data/%/genome.fa.gz, $(SPECIES))
all : $(GENOMES)

I have deleted Make equivalent since was not really related to the problem.

authors: طارق بك محمود


text:
How can I specify to DEseq2 to only perform comparisons on which I am interested with?
I am currently performing a large RNA-seq analysis from mice PBMCs. The dataset contains around 6,000 transcriptomic profiles and I would like to use DESeq2 to identify the sets of differentially expressed genes in the different conditions. In total, I have 100 biological stimulations, and for each stimulation I have 30 control samples and 30 samples treated with a molecule of interests (so is have: (30 controls+30 treated) * 100 stimulations = 6,000 samples).
I want to identify, for each stimulation, the sets of differentially expressed genes between control samples and treated samples. I do not want to compare samples from the different stimulations. In total, I would like to have thus 100 lists of differentially expressed genes.
I have started to use deseq2 to identify these lists but deseq2 is spending lot of time to perform comparisons on which I am not interested with (comparisons between the biological stimulations).
For now, I have a table sampleTable which looks like that:

I am using DEseq2 using the following command:
DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)

Could you please help me with that? How can I specify to deseq2 to not perform the comparisons between the biological stimulations, but rather between control and treated sample within each stimulation ?
Thank you and best,

I am not sure that you can do that... why are you not interested between control and treated of the different stimulations?

You can specify the exact comparisons you want in the results() function. So:
dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds = DESeq(dds)
res = results(dds, contrast("condition", "treatment1", "control1"))

Note that condition should not be, "stim001-control1, stim001-control2, etc.", but instead, "stim001control, stim001control, stim001control, stim001treatment, etc.". Don't put minus signs in your levels and don't number your samples in them (they're not in different groups.
The last command would be repeated for each of the comparisons. Note that this will still be quite slow due to the number of samples and the number of groups. With so many samples, you might just make a separate sampleTable for each of the comparisons you want to make. Each would then only contain the samples relevant for that comparison, so 60 total in each.
Alternatively, if you want to fit all of the groups at once then it's likely that limma/voom will prove to have better performance. It uses different math that happens to be quicker with large models.
BTW, if the slowness is coming from DESeqDataSetFromHTSeqCount() then you can mere everything on the command line into a single matrix. A couple lines of python or even the join command can do that quickly enough.


text:
How to proceed after contigs.fa file obtained from VELVET assembly
I want to assemble a linear genome and I have a reference genome too.
I ran the following commands for de novo velvet assembly:
velveth auto 31,45,2 -fastq -shortPaired1 -separate 2-CA-SVA-2_S2_L001_R1_001.fastq 2-CA-SVA-2_S2_L001_R2_001.fastq

I chose auto_33 for further analysis
velvetg auto_33 -exp_cov auto

Now my output files inside auto_33 folder are:
contigs.fa
Graph2
LastGraph
Log
PreGraph
Roadmaps
Sequences
stats.txt

But I have no clue how to move ahead. Can I do anything else with velvet or do I have to switch to another program?

Blat the contigs

Would you mind elaborating ?

How to proceed depends entirely on  what you want the assembly for. Are you looking for variants? Is this the genome of a related species which has not been annotated and you want to find its genes? Are you looking for one gene in particular? For syntenic regions? We can't tell you how to proceed if you don't tell us what you need. Please [edit] your question and clarify.

The contigs.fa file is a FASTA formatted file of the final contigs generated by velvet.  I would start assessing the quality of the assembly that was generated. Here are some questions to get you started.  Really the next steps depend on how you plan to use the assembly.
What is the distribution of contig lengths? Average and N50 contig length
How does the assembly compare to the reference?

How much of the reference does the assembly cover?

When you align the original reads back to the contigs how well are they covered?

Do any of the read alignments not make sense?
Are there sections of the contigs that are not covered?


text:
Using pythoncyc to get ec-code of a reaction
I am trying to collect the ec-code (and gene) of a reaction. I have tried the following:
import pythoncyc
pgdb = pythoncyc.PGDB("sco")
r = "RXN1A0-6128"
gene = "SCP1.233B"
enzyme = pgdb.enzymes_of_reaction(r)
enzyme = pgdb.enzymes_of_gene(gene)
pgdb.full_enzyme_name(enzyme)

However, I only get an enzyme id, and not the ec-code, which is visible on biocyc.

Have you looked at the documentation for pythoncyc?  I have not used that library before but it seems like the .full_enzyme_name() function would only return that.  Is there another function to get the ec-code?

It seems that ec_number is a key in the dictionary you get when you call PFrame() with getFrameData=True.

BioGeek is correct, however I will post a full answer.
import pythoncyc
# To get ec-code
db = pythoncyc.select_organism('sco')
reaction = db["RXN1A0-6128"]
ec_code = reaction["ec_number"]

# To get genes
pgdb = pythoncyc.PGDB("sco")
genes = pgdb.genes_of_reaction("RXN1A0-6128")


text:
Download data from the Human Microbiome Project via ascp
I have asked this question in biostars, but I am trying here as well as people working with "omics" data might be able to help. I think my issue relates understanding how large data storage on online server works. I am trying to download large data from the Human Microbiome Project. These are available on online server and can be downloaded via aspera (ascp).
I can't figure out the right path to get the data. On their website, they give an example of command to retreive an example dataset:
 ascp -l 100M -QT user@aspera.hmpdacc.org:/ReferenceGenomes/test_genome.tar /path/to/local/directory

I executed this, using my own username, from this I am asked my password, and I can authenticate. But then I get the following error:
(Error: Server aborted session: No such file or directory)

Then, I downloaded the manifest file via the HMP portal. In this manifest, they provide the path to which these files are supposed to be available. I tried to execute again the command with one of those path, provided by the manifest file:
ascp -l 100M -QT user@aspera.hmpdacc.org:/files/5a950f27980b5d93e4c16da124000378  /path/to/local/directory

I again get the same error of "No such file or directory".
The full url provided by the manifest
I tried putting different chunks of this address into to ascp command but nothing works. I think there is something I don't understand about how to access the server via aspera and what is expected to go in this path. When I copy paste this url in google, this directly download the file I want. So I'm close, but I need to do it via command line as I need to download hundreds of these files and they are huge (6-10GB each).
How can I download the data of the human microbiome project using ascp?

Is the manifest file a list of links like you provided? I mean to ask is the manifest file hundreds of lines with a URLs like http://downloads.hmpdacc.org/data2/hhs/genome/microbiome/wgs/analysis/hmwgsqc/v2/SRS013951.tar.bz2per line? It would helpful to see the format of the manifest file. The first 5 lines should suffice.

I can't test for ascp because they do not longer support my current setup but you can follow the FTP instructions. With FileZilla:
url: public-ftp.ihmpdcc.org

It looks like @CAsimonet has already gone through the instructions on the website and is have fine-tuned issues with getting a lot of specific files.


text:
5'UTR and 3'UTR annotation in yeast
I am working on a project in which I need to compute several parameters (such GC content and length) of 5'UTR and 3'UTR sequences of Saccharomyces cerevisiae yeast genes.
The problem is finding a proper annotation for these regions in yeast. I have tried with BiomaRt and these sequences are not available for yeast. I have also tried with UCSC table browser to obtain a BED file with 5'UTR or 3'UTR sequences but it is returning only coordinates of tRNAs, pseudogenes, ncRNAs.
In summary, is there a straight-forward way to obtain an annotation for 5'UTR and 3'UTR regions in Saccharomyces cerevisiae yeast?

Sorry, is for Saccharomyces cerevisiae. I am going to edit the question just in case.

isn't there an equivalent of AceView for Yeast ? https://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/index.html

I am unaware of any "official" or gold-standard UTR annotations in S. cerevisiae.
One option is to use the annotations from the TIF-Seq publication (Pelechano et al. 2013).
The GSE39128_tsedall.txt.gz file contains the major isoforms identified. It would be up to you to computational associate each transcript with a given gene. It is also up to you to determine which major isoforms are present in your dataset. Furthermore based on my own experience there is often a significant number of discrepancies between what I see in an NGS dataset (e.g. RNA-Seq) and what is called as a possible isoform in these annotations. However, a lot of this could be related to the inherent technological limitations of the method.
If you look at the SGD genome browser they provide some additional UTR annotations. Derived from Nagalakshmi et al. 2008.
You can find the relevant annotations for those UTR's in their supplemental data section. The file is 1158441_tables_s2_to_s6.zip and inside you will find TableS4 which seems to contain what you seek.
There might be more resources but these are the two main ones I'm aware of.

I found the following two files in https://downloads.yeastgenome.org/sequence/S288C_reference/:

SGD_all_ORFs_3prime_UTRs.fsa
SGD_all_ORFs_5prime_UTRs.fsa

According to the README files in the same directory, these are (the README for the 5' file is equivalent):

SGD_all_ORFs_3prime_UTRs.README
Information about the SGD_all_ORFs_3prime_UTRs.fsa file.
This file contains sequences for predicted 3' UTR sequences for all
  characterized, uncharacterized and dubious ORFs in S288c.
3' UTR coordinates were retrieved from YeastMine
  (yeastmine.yeastgenome.org) on July 3, 2017 using the stop coordinates
  of both the predicted transcript and gene sequences. A BED file with
  these coordinates was made. BED files were fixed on Sept 22, 2017 to
  account for strandedness. Sequences were batch downloaded again in
  FASTA format from USCS Genome Table Browser on Sept 22, 2017
  (https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=604915183_aUc9FZ1vwyAccGti1uI0wteYEM7x)
File created: 2017-08-25

Those therefore seem to be the 3' and 5' UTR sequences of all yeast genes in FASTA format and should be exactly what you need.

The link is broken

@Pengyao yes, apparently yeastdownloads.org are using a bad SSL certificate. I guess they'll fix it eventually.


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.


text:
Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)
I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.
I have a population allele count data, and it looks like this:
1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0

Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").
In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.
I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:
1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0

As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?
Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...
awk '{for(i = 1; i <= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i<=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i <= NF; i++)
        allele2[i] = $i;
            for(i = 1; i <= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i <= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'

It is strange that with R it takes so much time. Usually with matrices it is quite fast. Maybe you can post the code in R to translate it to python or perl, or improve it. perl is usually better for text editing, as with numbers it is a bit more difficult but it is quite fast usually. What have you tried in those languages? Did you translate your R code? Did you start from zero?

File is 8GB. I included bash solution but do not know exactly how to translate it into perl or python!

If you have a solution that works in on language then I dont' think this is the right forum.  It is no longer a bioinformatics question it is a programming question that happens to use biological data.

It seems like you do 7 loops in your awk code. I am not an expert of awk, but this calls for an optimization, I think it could be reduced to 4 or so. How is your code in R? (I am more confident of my R skills, and perhaps I could help you) And as Bioathlete says this might suit better to code Reviw. Where they could help you speed and improve your code

Could you also share the R solution for for those of us that find awk one liners to be completely impenetrable?

@ LIopis, OK. what is code review? is that what you meant? is that a forum?

@Anna, code review is when someone looks your code to improve it. See this site in the stackexchange network.

Your "bash" code is an awk solution

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:
allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000),
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

Edit
I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution
I used readr for file IO:
library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

Here are the timings:
cat Benchmarks/R.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.

Intelligent use of vector reusing! I would recommend scale  though

Thanks. I nicked the idea from here

Thanks for the benchmarks. Maybe you can answer this question of mine, about the meaning of the benchmarking values given by psutil (via snakemake): https://stackoverflow.com/q/46813371/1878788

I'm new new to benchmarking with Snakmake and also couldn't find any useful documentation for it. I did my best to answer with the info I've found, but I guess it's a case of the blind leading the blind

In python, if your files can fit in RAM, you might try using pandas and cytoolz.
import pandas as pd
from cytoolz import interleave

def odd(num):
    return num % 2

def even(num):
    return 1 - (num % 2)

# Read in two separate dataframes
alleles_1 = pd.read_table("snps.txt", skiprows=odd, sep="\s+", header=None)
alleles_2 = pd.read_table("snps.txt", skiprows=even, sep="\s+", header=None)
sums = alleles_1 + alleles_2
# https://stackoverflow.com/a/45565745/1878788
pd.DataFrame(list(interleave([
    alleles_1.div(sums).fillna(0).values,
    alleles_2.div(sums).fillna(0).values]))).to_csv(
        "freqs.txt", sep="\t", index=False, header=False)

I don't know how efficient this is. On my workstation it runs in 1m27 on data obtained by taking 1 million times the example given in the question. See heathobrien's answer for benchmarks on larger data.

Interesting, though I actually used the example times 10 million. I added some proper benchmarking to my anwser

@heathobrien Indeed, I had misread your code. With 20 million snps, on my workstation where Firefox is eating one third of the 16GB RAM, my script fails on a MemoryError after 4 minutes.

Here is an implementation in pure python. Probably not the fastest, but you should not encounter any memory issues with this.
Note: this is python3, use as follows: cat snp_counts.tsv | python allele_frequencies.py
#!/usr/bin/env python
"""
allele_frequencies.py
"""
import sys

is_odd = True
for line in sys.stdin:
    parts = line.strip().split()
    if not parts:
        continue
    allele_counts = [ int(p) for p in parts ]
    if is_odd:
        allele1 = allele_counts
        is_odd = False
    else:
        allele2 = allele_counts
        allele_sum = [ sum([a1,a2]) for a1,a2 in zip(allele1,allele2) ]
        allele1_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele1, allele_sum) ]
        allele2_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele2, allele_sum) ]
        print(*allele1_freq, sep = '\t')
        print(*allele2_freq, sep = '\t')
        is_odd = True

thanks for the code. I tried your code but I get an error message:  File "./allele_frequencies.py", line 21
print(*allele1_freq, sep = '\t')
^
SyntaxError: invalid syntax.

@Anna1364 Are you sure you are using python3?

Here's my quick Perl attempt after misreading the question, which divides each row elements by the row total:
cat popCountData.txt | perl -lane '
  my $total=0;
  grep {$total += $_} @F; ## create per-row total
  foreach (@F) {printf("%0.2f ", ($total) ? $_ / $total : 0)} @F; ## divide elements by total
  print ""; ## print a line break (note: line will end with a space)
'
0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.20 0.00 0.00 0.00 0.00
0.00 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.25 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.13 0.13 0.00 0.00 0.13 0.07 0.00 0.00 0.00 0.00 0.13 0.27 0.00 0.00 0.00 0.13 0.00

The challenge here is that the processing is in chunks of two lines, which means a couple of things:

The previous line (at least for even-numbered lines) needs to be stored
The sum of counts (or whatever statistic is necessary) will have different processing on each line

Unfortunately, I can't quite see how to get from your example input to your example output. For example, how does a count of 2 and 4 turn into a frequency of $1/3$ and $2/3$?

Because in that colum there is a 2+4 = 6 So 2/6 = 1/3 and 4/6 = 2/3, The scale is by each two rows and column wise.

@gringer, thanks but it does not give the output that I want. I added bash solution.Yes, LIopis is right


text:
BiomaRt error: Error in martCheck
I am experiencing an error while using getSequence() in biomaRt.
The code is given below:
library("biomaRt")
ensemblp <- useMart(biomart="plants_mart",host="plants.ensembl.org")
osj = useDataset("osativa_eg_gene", mart = ensemblp)
genes <- read.csv("osj.csv")
osj_id <- genes$ensembl_transcript_id

os_genes <- getSequence(id = osj_id, type = "ensembl_transcript_id", seqType = "peptide", mart = ensemblp)

The error message I am getting is :
Error in martCheck(mart, c("ensembl", "ENSEMBL_MART_ENSEMBL")) :
This function only works when used with the ensembl BioMart.This
function only works when used with the ENSEMBL_MART_ENSEMBL BioMart.

Did I make any mistake in code? Please help me fix the error. I tried asking this in Bioconductor support, but no one responded yet.

It doesn't seem to be a mistake in the code, but simply what you try to do isn't working in that BioMart

Did you find and review this post on bioconductor?  Does it help you? https://support.bioconductor.org/p/74761/

getSequence has only been enabled for the main Ensembl (vertebrates) biomaRt. This is all at the Bioconductor end. You'll need to use a getBM instead, for example:
os_genes <- getBM(attributes =c("ensembl_transcript_id", "peptide"), filters = "ensembl_transcript_id", values = osj_id, mart = osj)

Biomart seems to be down most of the time so i think you can parse the ensemble ID and names from the annotation file which i did as i was having similar problem .
You can use this to make a file with ensemble ID and gene names for you use
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
            printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt


text:
Why do I obtain very strange results with DESeq2?
I am using DESeq2 to perform a differential expression analysis, but I obtained very strange results for some genes. For some genes, I have very high log2FoldChange with very low p-value as displayed in the following table.

When I look the normalized expression values for these genes.  I can see that these genes are expressed in only one single sample (condition1= samples from 1 to 7, condition2 = samples 8 to 13).

Could you please tell me how to avoid DESeq2 to detect these genes as differentially expressed ?
I am using to following commands:
dds                      <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds                      <- DESeq(dds)
res                      <- results(dds, contrast=c("condition",condition1,condition2))

What are the size factors for these samples (sizeFactors(dds))? I suspect that sample10 is just really wonky and throwing everything off, so you'll probably want to just exclude it.

I have check that but the factor is in the range of other samples... I have similar situations with other samples..

Have you made any of the diagnostic plots related to cook's cutoff distance? If so, post some. If not, make some (see the here for an example).

I'm guessing here, but I suspect this is at least in part due to the fact that all your samples in condition 1 are not expressed at all.
Formally you have an infinitely large fold-change between condition1 and condition2. Most tools solve this by adding a small pseudocount to unexpressed samples, but this will still result in a very high fold change. So far for the big fold-change.
If I had to guess at the reason why this is significant I would think in terms of dispersion estimation. This Bioconductor post has some discussion with the authors of DESeq2 on how this works.
Finally: calculating differential expression on these genes in these samples will never be really informative, as they really are only expressed in one sample. The way I see it you have two options:

Exclude the sample from your analysis. Straightforward, but does not solve this problem for other samples.
Filter your genes before calculating differential expression, i.e. remove all genes that are not expressed in a large number of samples. This approach is standard in for instance sleuth. A good starting point would be to remove genes that are not expressed in 75% of your samples.


text:
Can structural - equation modeling be used in Bioinformatics?
I am a bioinformatics student and I was learning the basics of structural equation modeling (SEM).
I was wondering, are there any applications of SEM on  bioinformatics? Are there any strong sides for use of structural equation models using data on bioinformatics?
I did a quick search and found almost nothing. While it seemed like a learning algorithm, I can't say I figured it out completely.

I am not sure I understand your question, could you [edit] to clarify it? What do you mean by "strong sides for bioinformatics"? What do you want to know about bioinformatics? Or are you "just" asking about the algorithm of SEM? BTW here is a paper using SEM.

I mean, I haven't come across the usage of SEM in bioinformatics related articles. Bu it is widely used in education, questionaries and behavioral studies. In summary, I was wondering whether this method suitable for genomics studies or not.

I have came across them in gene regulatory network analysis before. see here and here. However I haven't seen them in other contexts, but that's not to say they don't exist.

SEM specialises in identifying 'latent variables'.
SEM are extensively used for questionnaire data. In questionnaire data, latent variables are the 'psychological' factors behind given (human) behaviours.
Applying SEM to explicit molecular mechanisms doesn't make much sense unless there is a hidden component in a mechanism that has not been identified. However, even here the SEM will not tell you what the 'latent variable' represents.
Thus I don't believe there is a good application to bioinformatics, outside human phenotypes derived from questionnaire data, which might be mapped on biological variables, but I could be wrong.


text:
bigWig to bed for regions above/below threshold
I'd like to threshold a bigWig file on a value and get all regions above or below that value in bed format. Ideally without too many file conversions. Are there any tools to do that?

Welcome to Bioinformatics. What have you found? Where did you look up ? (To avoid looking at the same sites) What have you tried? I think that we have some other questions in this site about bed and bigwig. Aren't they helpful?

cross posted: https://www.biostars.org/p/282410/

You can do that with a bit of python:
#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.

Convert bigWig to WIG via bigWigToWig:
$ bigWigToWig signal.bw signal.wig

Convert WIG to BED via BEDOPS wig2bed:
$ wig2bed < signal.wig > signal.bed

Filter BED by its score column with awk:
$ awk '($5 < 1)' signal.bed > signal.filtered.bed

The test $5 < 1 filters elements from signal.bed that have a score value less than 1. Change this to the threshold test you would want to apply.


text:
Interpreting large z scores from microarray data
I have microarray z scores for gene set enriched (GSE) pathways that range approximately from -3 to 14. I understand that an absolute value of 1.96 has a p-value = .05 and significant for our purposes. My question is whether I can consider extreme z scores as more significant?

This is a statistical question. What a higher z value means is that it is more extreme (if you assume a normal distribution), thus a more extreme value is less likely to happen by chance. Which is translated to a (lower) p-value.
What it is significant and not is arbitrary. Prior to the experiment you must decide which is your threshold, how much chance you allow. Apparently you decided to select 0.05 as your threshold, which means that one of 20 results might be due to chance. But your extreme z score is not more significant, it is less expected (under the null hypothesis, those p-values might change if you change your null hypothesis, ie you think your genes expression don't follow a normal distribution) but as significant as your 0.0499999 p-values.


text:
What does "motif 2 bp or longer" mean?
I'm reading a preprint on a new somatic variant discovery tool, and the first sentence of the Discussion section has me confused.

Across the four datasets analyzed in this study, we discovered that the major source of disagreement between callers originates from somatic variants called within STRs, in particular if the motif is two base pairs or longer.

The "two base pairs or longer" bit is what confuses me. I interpreted "motif" as referring to the short tandem repeats (STRs), which are by definition 2 bp or longer. A motif of < 2 bp is...a homopolymer, not really a motif in the common sense.
Am I missing something here? Or is this simply an awkward word choice on the part of the authors?

I doubt you’re missing anything, I suspect they’re using confused wording.

What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.


text:
RapMap: reference transcriptome for simulated reads
I am very new to bioinformatics and trying to repeat the benchmark in the RapMap paper with an experimental tool working in a similar but different fashion.
In the paper (taken from their github) about RapMap, an alignment/mapping tool for rna-seq reads to a reference genome, the dataset used for their benchmarks is described like this:

To test the accuracy of different mapping and alignment tools in a
  scenario where we know the true origin of each read, we generated data
  using the Flux Simulator (Griebel et al., 2012). This synthetic
  dataset was generated for the human transcriptome from an annotation
  taken from the ENSEMBL (Cunningham et al., 2015) database consisting
  of 86090 transcripts corresponding to protein-coding genes. The
  dataset consists of  48 million 76 bp, paired-end reads. [...]
  When benchmarking these methods, reads were aligned directly to the transcriptome, rather than to the genome.

I understand that they generated reads using flux simulator, but I don't understand which transcriptome they mapped them against. I can't find any other information in the paper about which annotation they used either.
I am assuming one would need to generate an entire transcriptome in .fa format from the GTF and genome? Or did they generate a transcriptome from only the simulated reads and then mapped the same reads against it?
GENCODE offers GTFs and FASTAs for gene-encoding transcripts here, might that be the kind of thing they used? Would the reads simulated from these GTFs fit to these transcript FASTA files?

FYI, I took the liberty of asking Rob Patro to have someone from his team reply.

Thanks for your interest in RapMap.
At that time we were using flux simulator for simulating read sequence data. We used the genome and gtf file together as an input to flux.
I dug into the scripts and got hold of the gtf link ftp://ftp.ensembl.org/pub/release-80/gtf/homo_sapiens/Homo_sapiens.GRCh38.80.gtf.gz,
and the genome file link ftp://ftp.ensembl.org/pub/release-80/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.
We used rsem-prepare-reference to create the reference transcriptome file later in order to run RapMap and other transcriptome based quantification tools. I hope this helps.
PS: I am a coauthor of the RapMap paper.

If they generated reads from the transcriptome then the started with the FASTA version of the ENSEMBL transcriptome and then aligned the reads back to it.  For a project like this I don't think that you need the GTF files, since you don't need to know the location of the transcript on the larger genomic reference.  If you can't find an ENSEMBL transcriptome then you can use the GTF and the bedtools program to create a FASTA.

The read simulator they used, Flux, can only generate reads from GTF format annotation and a genome, not from transcriptomes in FASTA format. For the time being I'm trying to work with the current GENCODE releases, which offer transcriptomes in fasta format, corresponding to their GTF annotations, but I will also look into your suggestion of using bedtools.

authors: Андрей Буфетов


text:
Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern
What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?

Are you interested in READs that start with a particular sequence or ALIGNMENTS that start with a particular sequence, because they are not the same thing for 2 reasons:

The sam format records sequences in genome orientation. Thus if the sequence of the read starts AAAA, but the read maps to the '-' strand the BAM seq field will not start AAAA, but end TTTT

Mappers may clip or otherwise only map part of a read. When they do this, they might only record part of the sequence of the read in the SEQ field. Thus a line with the CIGAR 5H90M will not contain the first 5nt of the read seq

First option
My solution uses the following steps:

use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C )
use jjs (java scripting engine) and the htsjdk library to build a bufferof reads having the same name. If any read starts with a pattern, print the whole buffer. The htsjdk should be already contained in the picard.jar.
sort back to coordinate order using samtools.

All in one:
java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
jjs -cp /path/to/picard.jar  script.js |\
samtools sort -T TMP -o out.bam -

with script.js
var out=java.lang.System.out;
var File = Java.type("java.io.File");
var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
var iter=samReader.iterator();
var buffer=[];
for(;;)
        {
        var rec=null;
        if(iter.hasNext()) rec=iter.next();
        if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
            {
            var i=0;
            for(i=0;i< buffer.length;++i)
                {
                if(buffer[i].getReadString().match(/^AAAA/)) break;
                }
            if(i!=buffer.length)
                {
                for(i=0;i< buffer.length ;++i)
                    {
                    samWriter.addAlignment(buffer[i]);
                    }
                }
            if(rec==null) break;
            buffer=[];
            }
        buffer.push(rec);
        }

samReader.close();
samWriter.close();

A second option is to use samtools view and awk to export the read names:
 samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt

and to use Picard FilterSamReads and READ_LIST_FILE to filter the reads.
UPDATE: 3rd solution:
I've updated my tool samjdk: there is now a new option --pair to work with query-name-sorted BAM files.
The command would be:
java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
java -jar dist/samjdk.jar --samoutputformat BAM --pair \
     -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
samtools sort -T TMP -o out.bam -

Does this account for the fact that BAM stores the sequence of a read in genome orientation, rather than read orientation? So if the sequence of a read mapped to the plus strand starts with AAAA it will be recorded in the BAM as ending with TTTT.

@ia no but it would be easy to test this: (R.getReadNegativeStrandFlag()?R.getReadString().endsWith("TTTT"):R.getReadString().startsWith("AAAA"))

I got it to work now with samjdk as explained in 3rd solution. Thanks a lot!

For at least one in a pair matching the sequence, if the read names of pairs are identical, some grep magic could do the job:

Print the sequence at the beginning of the line and read the name in the second column
Find all lines starting with the sequence
Extract only read names

And then just grep all found reads and save into the BAM file:
SEQ=<my_sequence>
samtools view unfiltered.bam | \
   awk '{print $10 "\t" $1}' | \
   grep "^$SEQ" | \
   cut -f 2 > list_of_matching_names
samtools view unfiltered.bam | \
   grep -f list_of_matching_names | \
   samtools view -b > filtered.bam

Two small notes:

You might want to consider also reverse complementary sequences.
You also might want to keep header of BAM file (you can extract it at the beginning and add it in the end).

I am one of the authors of the UMI-tools pacakge, and although you probably could hack together a way to do this with our package, its not what its designed for and the solutions above are probably better.
The best way to do this would be to filter the fastq before mapping if thats a possibility, if it has to be the BAM, then here is a python implementation similar to Pierres java solution above. This one fiddles about with strands and cigar operations a bit more.
First sort the reads by name with samtools sort -n
then
import pysam
inbam = pysam.AlignmentFile("mybam.bam")
outbam = pysam.AlignmentFile("mybam.filtered.bam", "wb", template=inbam)
seq_to_find = "AAGG"
rc_seq_to_find = "CCTT"

def test_and_out_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find):
    '''This function takes a list of reads and check if any of them
    have the prespecific sequences at their start. If the sequence
    is detected the whole buffer is output. The start of the read
    is the end of the alignment if the read is on the '-' strand.
    Also checks that the read doesn't start with hard-clipped sequence'''

    do_output=False

    for segment in read_buffer:

        if segment.is_reverse:
            # the last operation is stored in the first part of the
            #last of the cigar tuples. Operation 5 is Hard-clip
            output_condition = (not segment.cigartuples[-1][0] == 5) and \
              segment.query_sequence.endswith(rc_seq_to_find)
        else:
            output_condition = (not segment.cigartuples[0][0] == 5) and \
              segment.query_sequence.startswith(seq_to_find)

        if output_condition:
            do_output=True
            break

    if do_output:
       for segment in read_buffer:
           outbam.write(segment)

last_name=""
read_buffer=[]

# Loop over reads in the file
for read in inbam.fetch(until_eof=True):

    # if the name hasn't changed, the name is part of the
    # same "pair"
    if read.query_name==last_name:
        buffer.append(read)
        continue

    # if we got here then the a new pair is being examined
    # test the old pair and output if neccesary.
    test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)

    # make a new buffer with just the single new read.
    last_name=read.query_name
    read_buffer=[read]

# the last set of reds will not have been output.
test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)


text:
Filtering out all seqs with mutations of list2 from list1
I have 2 lists, list1 and list2, of protein sequences of the same given gene in different strains. In list1 are reported amino-acid sequences from resistant strains, and in list2 are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence.
I want to filter out (in Python 3) from list1 all of the variants that are also present in list2. The final aim is to have a list of the sequences that have unique mutations that are not in the 2nd list. The mutations can be of any type (e.g., indels, mismatches).
All the sequences are of at least 90% identity to some reference sequence  (filtered out the low matches using seq alignment).
Any enlightenment / package / tool for this job will be appreciated!
Toy example:
a, b, c : amino-acids
For this input:

Reference seq: aaabbb
Resistant strains: list1 = [aacbbc, ababbb]
Susceptible strains: list2 = [aacbbb, aaabbc]

The expected output would be: [ababbb], since the mutations in list2 cancels out the item aacbbc.
Of course this is a simple example and it is possible to have more complex scenarios (e.g., longer indels and mismatches).

Is there only one mutation per sequence entry in the lists?  Could you align the two lists to each other and look for perfect matches?

For the first step, I indeed filtered out the exact matches. But I also need to filter out other cases. For example if I have a sequence in list1 with 2 mutations and I have 2 sequences in list2 that each one of them has one of those mutations, then I want to filter this sequence from list 1 as well.

In which format is each list? It seems in a Fasta format but I am unsure because usually the mutations are not stored in these format.

The lists contain strings (extracted from transcriptome fasta files) of the amino acids that got score of at least 90% of the max score against some reference sequence.

I think the best thing to do would be to isolate the mutations themselves into a VCF and then there are tons of tools that do that filtering step.  Otherwise I think you will have to write something custom

Is there a python package that create VCF from 2 sequences? It's the first time I hear about this format and it sounds like the right direction.

Wait, so you have two lists of strings (sequence) and just want to find those strings present in file1 that are not in file2? What format are your sequences in? Fasta? One per line? Single letter amino acid code? Please show us an example of each of your two files and the output you would expect from that example.

All of those questions have been answered in comments.

added an example.

A VCF is usually output from a variant caller.  I would suggest aligning all the sequences to a single reference and then you can use samtools to report the differences in a VCF format.

@Bioathlete is there a python package that creates VCF? I use Biopython for alignments.

You will have to look at the alignment tool you are uses and see if it supports it directly or if it will write out a BAM file then you can use samtools to generate the VCF.

The point is that I want it to be a completely automated tool. As I understand it, samtools is an external tool and does not have a python API.

samtools is external but has a python interface via the pysam library.  I am suggesting that you either write something custom or you will have to branch out and use something in addition to BioPython.  The whole process can be automated via a single python script and even as a worst case you can use the python subprocess to access third party tools.

Ok, thanks. Last question: How can I generate a BAM file from 2 sequences in python? Couldn't find it in BioPython or any other package. Is there a package for that?

It would be an output format from the aligner you are using if it is supported.

Could you please give a name of package that has this functionality?

Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:
from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)

How can I get a list of differences between 2 proteins in a way that is convenient to work with? Is there a tool for that? It's not as simple as comparing strings (even after alignment). There are many possibilities (gap in one/gap in other/mismatch of different types/...)

Actually I think insertion deletion and mismatch are the only possibilities. Just store the position of the mutation in the reference sequence (or positions in the case of a deletion) and the sequence of the reference at that position/positions. You could also make a multiple sequence alignment of all sequences and discard sequences that from list1 that share non-reference characters with sequences in list2 without needing to keep track of reference positions, but that would require an external alignment program