text:
Removing genes with less than a correlation cut-off between two matrices
I have two matrices like this:
> head(DEG_log[,1:4])
            cl1       cl1       cl1       cl1
ALB    4.653796  7.046317  7.732135  7.758708
AQP9   7.798208  8.231112  8.734799  8.054087
CALML5 5.791524  9.058341  6.977369  9.503477
CCL4   9.739961  7.663303  8.934939  7.324958
CDK6   9.832115  9.710043 11.439345 10.330253
CHGA   5.139340 12.024606  5.988773  8.696836
>

I have the same patients and genes in two separate matrices coming from two sequencing panels with the same chemistry. Each sample has a match in another panel. I want to remove genes that have low correlation between two matched samples.
I have normalised and log-transformed the data by DESeq2.
I just plotted the correlation between each pairs of samples:
library("ggpubr")
ggscatter(my_data, x = "sample1", y = "sample1_1",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "onco", ylab = "biomarker)")

Like this:

by:
corr_test(a, b, 0.95)

I can calculate correlation, but I don't know how to remove genes with low or non-significant correlation between matched samples.

Welcome to the site! What have you tried? Did you found some errors or problems? Please along your wish describe your problem and what have you done so far to solve it, this way it is easier for us to suggest approaches you haven't tried or improvements on them. Do you have the same samples in both matrices (and ordered the same way)? Are the genes ordered the same way in both matrices ?

Thanks a lot, yes samples and genes are the same coming from same patients sequenced with same chemistry but separately

@FereshTeh is Mahta Mira another of your accounts?

Yes, Sorry. I am logging with Facebook, last night on my laptop this account did not work so I used that but now I am with office's computer

@DevonRyan Perhaps these accounts should be merged? Feresh(?) if you want it follow this guidelines. (Also tag me with @ llrs (without space) if you want me to receive the notifications)

@llrs Yup, I'll look into merging the accounts, it'll be a first for me.

Nice data and good question.
It is called "outlier" analysis or the analysis of the residual. Its an easy analysis, but . You calculated the standard deviation using the regression slope as the 'mean' and retain data within 1.96 standard deviations. You can plot the data like this and if you are lucky will see the residual forming a nice 'normal distribution'. Outliers outside 1.96 sd can be formally excluded.
Calling (and plotting) the residual from the regression slope will be possible in R, but I use external packages because its easy. BTW you should also try Spearman's coefficient to check you are getting the same answer.

Thanks a lot actually I need my data cleaned from genes with low correlation  because I want to combine my data sets. I was not able to find a line of code to remove such genes although I can plot them in correlation plot.

The specific statistic is "standardisation", its easy to calculation, but really its just working out the SD. Any could good stats package will  calculate residual, I like "DataGraph" (Mac).


text:
SNP vs Point Mutation
What is the difference between a Single Nucleotide Polymorphism (SNP) and a point mutation? I am quite confused in understanding these term as both of them refer to one base difference from the reference sequence?

The difference between the two depends on to whom you talk ;)
You are right: both refer to one base difference from the sequence. Usually the term "mutation" is used if the change has an impact on the phenotype. The "P" in "SNP" means "polymorphism". A lot of people use this term if more than 1% of the individual in a given population have this variant.
My recommendation: only use "mutation" if you are sure that your variant has an important impact on the phenotype. Only use "SNP" if you are sure that it is a commonly occurring variant. In any other case just use "variant". It is slowly becoming more common to talk about "SNV" (Single Nucleotide Variant) than about "SNP".

I would add that "mutation" is best used when it has an important negative impact on the phenotype. There's a lot of baggage that goes along with the word "mutation", so even when we talk about patients we usually prefer "deleterious variant" or something like that.

Personally, I only use "mutation" when talking about it actually happening. So somatic mutations are mutations, or de novo mutations where the child had a variant but the parent doesn't are mutations. But if it's germline, I always call it a variant. The other advantage of saying "variant" over "SNP" is that "variant" can be used to describe insertions and deletions larger than one base.

@DevonRyan But consider "silent mutation" or "adaptive mutation", these are not negative in their impact and people studying molecular evolution would not consider that the default; indeed most mutations are neutral

I think SNP is only really used at the population level while a point mutation may imply either population level (SNP) or between-species (substitution). Also sometimes, a point mutation is used for a single base pair insertion/deletion (although personally I wouldn't use it like that), while SNP is never an indel

@Chris_Rands I suspect this comes more from my human genetics background, where we worry about how patients will interpret what we're saying.

@Chris_Rands When I was working on comparative genomics and sequence evolution, I and those around me would use "SNP" for any form of single nucleotide change. Using the term specifically for common ones is something I only saw when I started looking at NGS data in a clinical context. I think this is one of those things that really change depending on your specific sub-field.

source: Stack Exchange

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

sort: votes


text:
How do I resolve disagreements in the determination of significant genes?
I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.
This is my data, log transformed RNA-seq:
              TRG    CDK6 EGFR  KIF2C CDC20
    Sample 1  TRG12  11.39 10.62  9.75 10.34
    Sample 2  TRG12  10.16  8.63  8.68  9.08
    Sample 3  TRG12   9.29 10.24  9.89 10.11
    Sample 4  TRG45  11.53  9.22  9.35  9.13
    Sample 5  TRG45   8.35 10.62 10.25 10.01
    Sample 6  TRG45  11.71 10.43  8.87  9.44
    ...

23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).
1- I tested each of 23 genes individually in this code and each of them gives p-value < 0.05 remained as a good predictor; For example for CDK6 I have done
> glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
>
> summary(glm)

Call:
glm(formula = TRG ~ CDK6, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.8286  -0.9317   0.6382   0.9500   1.5966

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  26.8950     9.6627   2.783  0.00538 **
CDK6         -2.6932     0.9765  -2.758  0.00582 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 64.062  on 54  degrees of freedom
AIC: 68.062

Number of Fisher Scoring iterations: 5

Finally I obtained five genes and I put them in this model:
final <- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))

but
> summary(final)

Call:
glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.5696  -0.9090   0.2233   0.7125   2.0131

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   3.2258    13.3877   0.241   0.8096
CDK6         -2.6101     1.0860  -2.403   0.0162 *
CXCL8         0.8702     0.9783   0.889   0.3737
IL6           0.9591     0.8883   1.080   0.2803
ISG15         0.9574     0.5101   1.877   0.0605 .
PTGS2        -0.4623     1.1005  -0.420   0.6744
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 53.917  on 50  degrees of freedom
AIC: 65.917

Number of Fisher Scoring iterations: 6

I am just seeing CDK6 as significant. I done cross validation
> cv.glm(df, final, K=nrow(df))$delta
[1] 0.2130050 0.2125748

Likely the delta values are not greatly differ (Model is right ???)
Then I tried Sensitivity and Specificity
> threshold=0.5
> predicted_values<-ifelse(predict(final,type="response")>threshold,1,0)
> actual_values<-final$y
> conf_matrix<-table(predicted_values,actual_values)
> conf_matrix
                actual_values
predicted_values  0  1
               0 16  7
               1 10 23
> sensitivity(conf_matrix)
[1] 0.6153846
> specificity(conf_matrix)
[1] 0.7666667

2- I went this way for step-wise method
full.model <- glm(TRG ~., data = df, family = binomial)

library(MASS)
step.model <- full.model %>% stepAIC(trace = FALSE)
coef(step.model)

> coef(step.model)
(Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17
-23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403
     MAGEA1         OSM       PTGS2
   496.6933   -637.0583   -500.1040

In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(>|z|) > 0.05
> summary(step.model)

Call:
glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 +
    ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial,
    data = df)

Deviance Residuals:
       Min          1Q      Median          3Q         Max
-1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -23770.6  2922188.9  -0.008    0.994
CDK6           -469.4    64920.6  -0.007    0.994
CHGA            376.4    46191.2   0.008    0.993
CSF3           -238.6    30857.1  -0.008    0.994
CXCL8          1404.6   172558.7   0.008    0.994
ERBB2           590.3    71989.5   0.008    0.993
IL6             789.5   100130.7   0.008    0.994
ISG15           842.2   102934.0   0.008    0.993
KRT14           423.0    52565.5   0.008    0.994
KRT17          -651.3    79489.6  -0.008    0.993
MAGEA1          496.7    60319.6   0.008    0.993
OSM            -637.1    82121.8  -0.008    0.994
PTGS2          -500.1    73541.7  -0.007    0.995

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7.7347e+01  on 55  degrees of freedom
Residual deviance: 1.4741e-07  on 43  degrees of freedom
AIC: 26

Number of Fisher Scoring iterations: 25

3- I went to lasso regression
> library(glmnet)
> lassoModel <- glmnet(
+   x=data.matrix(df[,-1]),
+   y=df$TRG,
+   standardize=TRUE,
+   alpha=1.0,
+   family="multinomial")
> #Derive coefficients for each gene to each subtype
> co <- coef(lassoModel, s=idealLambda, exact=TRUE)
> co
$`TRG12`
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept)  4.441322326
ALB         -0.010867331
AQP9        -0.051908908
CALML5       0.126359364
CCL4        -0.059869284
CDK6         0.283781437
CHGA        -0.088283426
CREB3L3      0.109988758
CSF3        -0.012311176
CXCL5       -0.030532485
CXCL6       -0.063610112
CXCL8       -0.084385732
ERBB2       -0.102303095
FGFRL1       0.148759702
IL1B        -0.072225606
IL6         -0.122671026
ISG15       -0.297386587
KLK5         0.039891341
KRT14       -0.004536929
KRT17        0.041693665
MAGEA1      -0.154075458
OSM          0.028987648
PTGS2       -0.034960405
S100A7A     -0.008190105

$TRG45
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept) -4.441322326
ALB          0.010867331
AQP9         0.051908908
CALML5      -0.126359364
CCL4         0.059869284
CDK6        -0.283781437
CHGA         0.088283426
CREB3L3     -0.109988758
CSF3         0.012311176
CXCL5        0.030532485
CXCL6        0.063610112
CXCL8        0.084385732
ERBB2        0.102303095
FGFRL1      -0.148759702
IL1B         0.072225606
IL6          0.122671026
ISG15        0.297386587
KLK5        -0.039891341
KRT14        0.004536929
KRT17       -0.041693665
MAGEA1       0.154075458
OSM         -0.028987648
PTGS2        0.034960405
S100A7A      0.008190105

>

Likely none of genes is significant so I just tried 5 genes from first glm
> finalLasso <- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link="logit"))
> summary(finalLasso)

Call:
glm(formula = df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.6780  -0.6733   0.1643   0.6925   1.9182

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4178    14.7456  -0.028   0.9774
CDK6         -2.2672     1.1537  -1.965   0.0494 *
IL6           1.2524     0.9197   1.362   0.1733
ISG15         1.1763     0.5479   2.147   0.0318 *
CXCL8         0.5413     0.6443   0.840   0.4008
CALML5       -0.7890     0.4219  -1.870   0.0615 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 49.839  on 50  degrees of freedom
AIC: 61.839

Number of Fisher Scoring iterations: 6

>

That 2 genes are significant; really confusing :(
Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not

First confusion: Why when I am putting each of 23 genes individually
in glm model, 5 of them returns p-value < 0.05 but when summary of
final model containing these 5 genes says only p-value of one of them
is < 0.05
Second confusion: When I am using stepwise regression, summary of
model says none of my genes is significant while some of them had been
remained as predictors by coef(step.model)
Third confusion: CDK6 as my significant genes in two of model shows
negative coefficient while I know that show be correlated with
dependent variable positively

confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value < 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is < 0.05.
Ans. Pr(>|z|) simply suggests that whether z-value deviates beyond average z-value and if so,what is the extent of its probability. less than .05 indicates a 95% confidence in computed z value i.e. an acceptable significance level. The summary of 5 genes simply tells significance of over-all model.

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


text:
How is consensus alignment for OLC assembly usually implemented?
Background
In C++ I want to implement a simple class that performs overlap layout consensus assembly but I can not figure out the most logical data structure to use for the consensus alignment step.
This is the step where all of the reads from the overlap graph are aligned and the consensus is called (slide 48). An example alignment is below. Asterisks above the column indicate some disagreement among all reads. Alignment position is shown above and below for convenience.
                         *    * *         *
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012
read 1   TAGATTACACAGATTACTGAT TGA
read 2      ATTACACAGATTACTGACTTGATG CGTAAACTA
read 3                   tTGACTTcATGGCGTAA CTA
read 4                     GACTTGATGGCGTAA CTAAG
read 5                           ATGGCGTAA CTAAGATCG
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012

Question:
Q1: What is the most logical way to implement the consensus alignment step in a OLC assembler?
Q2: What are the pros and cons of this and other possibilities?
Ideas
Idea - Vector of vectors. One column per position.
My initial idea for a data structure to store the alignment is to use one vector per read to store the read's character at that position. We would end up with j arrays (number of reads) all of length i (the total length of the overlapped alignments). This, however, would be space-inefficient if the locus were large and many positions in the vectors would contain a null value. In the example alignment above the first 23 positions of read 5 would be null values.
Idea - Graph. One node per position.
In another implementation idea, each position in the alignment would be represented by a node in a directed graph. Each node would have a property keeping track of the nucleotide counts at each position and might keep the read id's stored as well if there is a good reason to do so.

Implementing data structures and algorithms from scratch is a great way to build understanding and intuition. But if your goal is to create production-ready software, you may consider API calls to a mature third-party library such as fermi-lite.

For short reads, mainstream assemblers don't have a consensus step. The contig path directly spells the consensus sequence.
For long reads, the most popular way to derive a consensus is via partial-order alignment (POA). This is used by error correctors (similar to consensus) like pbdagcon and falcon_sense, and consensus tools like canu, racon, wtdbg2 and nanopolish. Basically, you start with an empty directed acyclic graph (DAG), align each read in turn to the graph and augment the DAG by inserting unaligned segments. When you align all reads onto the DAG, the most weighted path is taken as the consensus. You should read Jared Simpson's blog on this topic. It is much more detailed than "basically". You should also read the papers on all the mentioned tools above if you want to implement by yourself.
An alternative approach is Gene Myers' realigner (PDF here). It is not often used, though.


text:
What is a sensible forcefield choice for membrane proteins when using PDB2PQR?
I am generating PQR files for a membrane protein that is almost entirely buried in the membrane. The goal is to calculate the electrostatic charge across the surface of the protein. There are no lipids in the structure.
PDB2PQR has predefined options of amber, charmm, parse, tyl06, peoepb and swanson. It also allows user supplied forcefields.
The protein is an ion channel in the plasma membrane of a human. I'll be processing in memAPBS. However, I would greatly appreciate a "rule of thumb" answer as I plan to use the same method on TMPs from other organisms in various subcellular membranes.
What is a sensible forcefield choice for membrane proteins when generating PQR files for electrostatic analysis?


text:
Calculating Z-score from logCPM values using edgeR
I have the raw counts for RNA-Seq data. I converted counts data to logCPM using edgeR package.
Lets say I have a dataframe A with 15000 genes as rows and 100 samples as columns with counts data. I have used following code from edgeR tutorial.
logCPM <- cpm(A, prior.count=2, log=TRUE)
head(logCPM)[1:5,1:5]
               Sample1          Sample2          Sample3         Sample4
A1BG          -1.019429        0.8333739        -1.192960       -0.6717287
A2M            7.622280        7.8057767         7.932570        8.6409077
A2ML1          1.004541        1.0909119         5.462613        6.4149348
A4GALT         2.126502        1.8403954         5.322601        2.9469345
AAAS           5.785888        5.1639501         5.259789        5.3248026
               Sample5
A1BG          0.3121975
A2M           9.1949911
A2ML1         9.4740594
A4GALT        2.9214979
AAAS          4.3708063

After converting counts to logCPM I see there is another command for scaling which is used for heatmaps.
logCPM <- t(scale(t(logCPM)))

It is given there that scaling each row (each gene) to have mean zero and standard deviation one in page 15 edgeR
Is this nothing but z-score? or Is there any otherway to calculate z-score from logCPM data?
thanq

Yes this is nothing but a z-score. The scale function calculates column z-scores, but since you use t (transpose), it is actually performed on rows. There are many other ways to calculate z-scores, but if done correctly they will all give the same results.


text:
Sort vcf by contig and position within contig
Due to tabix constraints, I need to sort a vcf so that contigs and then positions within contigs occur in numerical order in the vcf.  I don't know if the following will sort positions within contigs, but I tried this:
cat in.vcf | vcf-sort -c > out.vcf

but was informed:
Old version of sort command installed, please run without the -c option.

I've read that can cause problems here. Further, my contigs are't numeric, but are something like this: 'Erumb1_s00000030'.
I see two possible solutions.  (1) I am on Mac OS Mojave, and if there were a way to install the a new unix-style sort, I could use it. I tried homebrew and didn't find anything. (2) there are likely some bash commands that would do it.  If anyone has an idea in either of these directions, I would appreciate it.

It might be worth posting a question about these tabix constraints. We could be able to help you fix that which would allow you to use bcftools.

You could simply use bcftools sort for it:
$ bcftools sort input.vcf > output.vcf

If you really want to use bash only, you can do this:
$ grep "^#" input.vcf > output.vcf
$ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf

The first command will write the header information to the new vcf file.
The second will sort by contig name and position and append the result to the new vcf file.
fin swimmer

I get this: Writing to /tmp/bcftools-sort.7BfymN
[W::vcf_parse] Contig 'Erumb1_s00000030' is not defined in the header. (Quick workaround: index the file with tabix.).  So, this is where I started, having to index the .vcf.  Looks like the grep approach worked, but I've got to give it a closer look.

@PeterPearman first check that your sort implementation deals with the -k flag correctly. The command in this answer should work fine for GNU sort (the default on Linux) but I'm not sure it works correctly with BSD sort (default on macOS). The man page I found seems to suggest that the V and g modifiers aren't supported. However, you don't really need the V and you can use n instead of g: sort -k1,1 -k2,2n.

bcftools is very strict about the correct vcf structure. It heavily relies on the header data. This should include contig information and any key used in the INFO and FORMATcolumn. You should always have this information there. bcftools is a very powerful tool. So be strict to your vcf file, that you can use it.

@terdon is right. V isn't needed. I'm used to it, handling human chromosome names. Doing it that way chr2 is located before chr10. But in vcf file the only conditions are, that the contig names build blocks and inside this blocks the position is ascending.

Use a subshell. First print the header and then the variant lines, sorting by chromosome and position:
(grep ^"#" input.vcf; grep -v ^"#" input.vcf | sort -k1,1 -k2,2n) > sorted.vcf


text:
Non-random access on a fastq file
There was a nice question on random access on a fastq file, however I have the opposite problem. I want to programatically access non-random fastq entries quickly (using C++).
Say I have several pairs of fastq.gz files that could be hundreds of millions or > 1 billion read pairs total. I know the number of read pairs in each file and I have a list of 10-20 reads that I want to access from each file. I only know the 0-based index of the fastq entry that I want to access, not the read ID. So if I want the  234th and 863rd reads the naive solution is to open the fastq file and scan through 864 records. I need to perform some tasks with these reads, then I move on and access/process another batch of reads.

I need to perform the above task tens of thousands of times, so the options are:

Low RAM usage. Many file IO ops. No disk usage: Opening/scanning/closing the fastq.gz files then processing the desired reads tens of thousands of times.

Extreme RAM usage. One round of file reading. No disk usage: Read in each fastq.gz file once, store all of the desired reads in RAM and access and process those reads as needed

Low RAM usage. One round of file reading. Some disk usage: Coming up with an intermediate file format similar to a bam that will allow me to parse and sort the reads with their groups and other metadata later on.

If it makes a difference the fastq.gz files are on spinning disks, not SSDs.

Load all lists into RAM and sequentially scan all FASTQ files.

You don't have a RAM problem because you only have about:

10-20 reads that I want to access from each file.

Please go for (2). If you are concerned about RAM usage, you may selectively write some of your data back to disk.
Tens of thousands of 10-20 reads is not a problem for modern computing.


text:
group samples based on shared mutations in a single multi samples vcf file
I am learning about vcf file formatting and have a multi-sample (>300) vcf file and aim to group (take) samples with shared mutations.
Can someone suggest a tool/command/script solve this problem?

FYI: Cross posted on biostars (https://www.biostars.org/p/358896/)

@finswimmer Cross posting on biostars is alright.

@SmallChess I'm not against cross posting, but one should mention it, where they ask the question as well. Doing this every one can see if there is already a solution.

As I asked I am learning and asked for help :)

@finswimmer posted this response on Biostars $ bcftools query -i 'GT="0/1"' -f '%CHROM\t%POS\t0/1\t%REF\t%ALT\t[%SAMPLE,]\n' input.vcf|sed 's/,$//'|awk -v OFS="\t" '{counts = split($6, samples, ","); print $1, $2, $3, counts, $4, $5, $6}' > stats_01.csv

I think the question needs to be much more specific for that answer to be useful for anyone else. What does 'shared mutation mean'? What does 'group' them mean? I think a small sample of the input and desired output is needed.


text:
Prediction of prokaryotic origins of replication (ORI)
I want to predict origins of replication (ORI) on hundreds of prokaryotic genomes. The most straight-forward solution would be to use most commonly used tool, Ori-Finder.
It uses integrated gene prediction, analysis of base composition asymmetry, distribution of DnaA boxes, occurrence of genes frequently close to oriC regions and phylogenetic relationships to predict ORIs.
However, I came across a paper by Parikh et al which describes a single measure (correlated entropy measure / CEM) with which they achieved 100% ORI-identification on 500 genomes.

I like the CEM approach more than a close-source Windows-only black-box.
Parikh tested their algorithm on NCBI genomes. How are ORI annotated on these reference genomes? I couldn't find such an annotation. Does the sequence simply begin at the ORI in leading strand direction by convention? Were the ORI of the reference genomes identified by experiment?

Hi and welcome to the site! I removed your second question since that is not something we can give a definitive answer to and this site only deals with questions that can be answered without relying on opinion. Sorry about that!

I have not seen the papers but I know the technique SVM. This is a machine learning approach its pretty good, but has been around this the 1980s.
Its quite simple known ORIs are given with the associated genetic trait, e.g. GC skew, SVM will then regress the GC skew of a known ORI against the GC content of the genome, there will be a training, parameter optimization and test split in the data, typically 80:20:20. De novo data will then be feed into the trained algorithm, in the same way the test data was, and a prediction will be generated. There are weaknesses, the authors have used only one ML technique, typically numerous will be used e.g. SVM, random forest, lasso regression etc... In addition, the prediction accuracy for GC is low (75%), in addition 500 genomes isn't great for ML, you need 1000s.


text:
Circos plot with circlize
This is the post where they made circos using R I'm just following that the part answered by Zuguang Gu, not not sure what I'm doing wrong my file, and the code:
library(ComplexHeatmap)
library(circlize)
 d11 <- read.csv("DC3_7_MUT.txt",stringsAsFactors = FALSE)
a <- as.data.frame(d11)
head(a)

ex <- melt(a, id.vars=c("Gene"))
head(ex)
dim(ex)

bb1 <- ex %>% na.omit() %>%
  group_by(Gene,value) %>%
  mutate(new = n()) %>%
  data.frame()

#head(bbb)
#dim(bbb)

bb1[bb1==""]<-NA

df <- na.omit(bb1)

bbb <- df

features1 = unique(bbb[[3]])
features_col1 = structure(rand_color(length(features1)), names = features1)

patients_col1 = structure(1:187, names = unique(bbb[[2]])) # I don't understand this

col1 = ifelse(bbb[[4]], features_col1[bbb[[3]]], "#FFFFFF00")
col1 = gsub("FF$", "80", col1) # half transparent
features_count1 = tapply(bbb[[4]], bbb[[3]], sum)
head(bbb)
# set color to white if it only exists in one patient/cell
col1[features_count1[bbb[[3]]] == 1] = "#FFFFFF00"
pdf("CIRCOS.pdf",20,25)
chordDiagram(bbb, col = col1, grid.col = c(features_col1, patients_col1))

dev.off()

The error

chordDiagram(bbb, col = col1, grid.col = c(features_col1,
  patients_col1)) Error in abs(df$value1[i]) :    non-numeric argument
  to mathematical function

I'm doing something wrong but I'm not able to find out, secondly I didn't get the idea about the usage of structure the one I commented.
Am I using the right dimension for the structure argument?
A bit of progress with a bit of tweak in the code
d11 <- read.csv("DC3_7_MUT.txt",stringsAsFactors = FALSE)

blank2na <- function(x){
  z <- gsub("\\s+", "", x)  #make sure it's "" and not " " etc
  x[z==""] <- NA
  return(x)
}

a1 <- data.frame(sapply(d11,  blank2na))

df <- a1
head(df)

df2 <- df[, colSums(is.na(df)) != nrow(df)]

head(df2)
dim(df2)

a <- as.data.frame(df2)
head(a)
ex <- melt(a, id.vars=c("Gene"))
head(ex)
dim(ex)
ex <- with(ex, table(paste(Gene,variable,sep='|'), value))
ex = as.data.frame(ex, stringsAsFactors = FALSE)
head(ex)
ex <- ex[which(ex$value != ""),]
head(ex)
dim(ex)
ex <- (ex[which(ex$Freq == 1),])
head(ex)
dim(ex)

features = unique(ex[[2]])
features_col = structure(rand_color(length(features)), names = features)
patients_col = structure(1:280, names = unique(ex[[2]]))
col = ifelse(ex[[3]], features_col[ex[[2]]], "#FFFFFF00")
col = gsub("FF$", "80", col) # half transparent
features_count = tapply(ex[[3]], ex[[2]], sum)
# set color to white if it only exists in one patient/cell
col[features_count[ex[[2]]] == 1] = "red"

#And the final chord diagram:

pdf("CIRCOS.pdf",20,25)

#png("circos.png",width =20,height=20,unit="cm",res=1000)
#chordDiagram(ex, col = col, grid.col = c(features_col, patients_col))

chordDiagram(ex, col = col, grid.col = patients_col, directional = 1,
             transparency = 0.25 , xmax = NULL,
             symmetric = FALSE, keep.diagonal = FALSE,
             direction.type = "diffHeight", diffHeight = convert_height(2, "mm"),
             reduce = 1e-5, self.link = 2,big.gap = 10,
             preAllocateTracks = NULL,
             annotationTrack = c("name", "grid", "axis"),
             annotationTrackHeight = convert_height(c(3, 2), "mm"),
             link.border = NA, link.lwd = par("lwd"), link.lty = par("lty"),
             link.sort = FALSE, link.decreasing = TRUE,scale = TRUE,
             link.arr.lwd = par("lwd"), link.arr.col = par("col"),
             link.largest.ontop = FALSE, link.visible = TRUE)

dev.off()

What exactly are you trying to achieve? What does each connection in the circos plot represent?

this i finally discarded ...


text:
How to set the position of groups in a Seurat object on a FeatureHeatmap plot
I am analysing singe cell sequence data and I have followed this tutorial, https://satijalab.org/seurat/pbmc3k_tutorial.html to perform QC and various differential analyses using the Seurat package on my data up till now.
How can I rearrange on the resulting plot, the order in which the groups of my samples appear when I use the FeatureHeatmap() function to visualise the differential gene expression of Gene1 and Gene2 across the cell clusters?
FeatureHeatmap(object = object, features.plot = c("Gene1", "Gene2"),
               group.by = "Conditions",
               pt.size = 0.25, key.position = "top",
               max.exp = 3)

produces:

but I want:

I don't think this is possible in Seurat v2, but in v3 you can change the factor levels of the grouping variable to change the plot order:
library(Seurat)

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

Change the order:
pbmc_small$groups <- factor(pbmc_small$groups, levels = c('g2', 'g1'))

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

Thanks a lot. However the function VariableFeatures is not present in my Seurat package. I get the error: Error in VariableFeatures(object) :    could not find function "VariableFeatures" @TimStuart

Yes, that function is part of Seurat v3. This answer will only work for v3

How do I get the Seurat V3 package? My computer is having issue downloading it automaticaally, so I wish to download the Seurat V3 package manually and then load it. Thanks. @TimStuart

Will it also be advisable that I move to Seurat V3 at this level or I may need to restart my analyses afresh and repeat all steps with Seurat V3?

You can install v3 using devtools: devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0'). I don't know how much effort it is for you to repeat analyses with Seurat v3, but you can also just update the Seurat object using UpdateSeuratObject. See here for more information: https://satijalab.org/seurat/essential_commands.html

I obtain this: devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0') : package ‘curl’ successfully unpacked and MD5 sums checked: Error: (converted from warning) cannot remove prior installation of package ‘curl’. UpdateSeuratObject(object): Object representation is consistent with the most current Seurat version.
The FeatureHeatmap() actually uses the alphabetical order of the conditions which I provided in the metadata slot. Is there a possibility to reorder the factor levels without harming the integrity of the data in the Seurat object?

Can you confirm Seurat v3 has actually been installed and loaded before running UpdateSeuratObject?


text:
Extracting some part of a list
I have two normalized gene expression values (log2 of cpm)
> head(biomarker[,1:5])
               A2        A3        A4        A6        A7
A2M    12.4071618 12.631601 12.889748 11.842511 11.574134
ABCB11  0.5151957  6.176788  5.414905  7.134915  1.247590
ABCC2   6.7244303  6.487794 10.178615  6.384132  6.089462
ABCC6   6.3977122  5.823627  7.370761  7.397091  6.071587
ABCF1   9.1609847  9.893258 10.116638 10.225520 10.838486
ABCG2   7.6841874  7.738293  6.833152  6.888041  6.030032
>

And
> head(immune[,1:5])
              A2        A3        A4        A6        A7
A2M    10.904748 11.388404  9.910614 11.513439 12.963609
ABCB11  5.011380  6.359443  7.145992  8.451947  7.722605
ABCC2   5.040461  6.477014  3.873996  6.409777  9.133971
ABCC6   7.798441  7.601848  9.948072 11.628533 12.701460
ABCF1   8.553597  8.615425 11.145903 10.289098 11.444140
ABCG2   6.224294  5.629375  8.293416  7.979859  8.603793
>

> dim(biomarker)
[1] 719  56
> dim(immune)
[1] 719  56
>

I have matched samples in these experiments and same gene too. I could roughly say in both experiments everything has been constant. However, I want to do Student's t-Test to find genes that they had inconsistent gene expression patterns between the two data sets.
By this function I am getting a list but I don't know how to extract the names of genes with p-value > 0.05 or whatever cut-off
> f <- function(x,y){
+     test <- t.test(x,y, paired=TRUE)
+     out <- data.frame(stat = test$statistic,
+                       df   = test$parameter,
+                       pval = test$p.value,
+                       conl = test$conf.int[1],
+                       conh = test$conf.int[2]
+     )
+     return(out)
+ }
> t=sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Any help please?
  str(t)

    [list output truncated]
     - attr(*, "dim")= int [1:2] 5 719
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:5] "stat" "df" "pval" "conl" ...
  ..$ : NULL
    >

> sapply(t, "[[", "pval")
Error in FUN(X[[i]], ...) : subscript out of bounds
>

> pvalueGenes <- sapply(colnames(biomarker),function(i) t.test(biomarker[i, ],immune[i, ], paired = TRUE)$p.value)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mean.default(x) : argument is not numeric or logical: returning NA
2: In if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :

 Show Traceback

 Rerun with Debug
 Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed

Did you try iterate by colnames instead? Have you searched the internet about how to get names of the variables you iterate over them? BTW you might be interested in the broom package, function tidy.

Sorry I don't know how to do that

A simple google search would give you answers. You could use tibble::rownames_to_column(biomarker, col_name). You could even do biomarker$col_name <- colnames(biomarker)

Sorry what is col_name please?

@FereshTeh you don't need to apologise for not knowing something. "col_name" is an abbreviation for column name. Use the help pages of those functions Ram mentioned to understand what do you need and what do they do

@FereshTeh Did you read and understood the help page of the rownames_to_column function?

Actually, no; You know, likely I have done t-test to find differentially expressed genes. I have 719 genes and two sets of samples while the function returns a list with 56 * 4 dimension that means t-test has been done between samples while I need to extract genes with p-value > 0.05 among two sets of datasets, that is why I got confused and I don't know rownames_to_colnames would solve which of these

@FereshTeh Ok, that is something you should [edit] into the question.

If you're going to downvote, add a reason why. Don't be passive-aggressive. It's bad for the community. Thanks!

@AlexReynolds It is not required to disclose, but I think that the downvote button already explain this: "does not show research effort, it is unclear or not useful". Asking for an explanation from downvoters is also passive-aggressive...

@FereshTeh I used col_name as a name for the new column that will be generated by the function rownames_to_column. I should have used 'col_name' (with quotes), sorry about that!

@AlexReynolds I kinda agree with llrs. In addition to not showing research effort, OP has a habit of editing her posts extensively and changing the entire premise once she receives a response, and such behavior needs to be discouraged.

Starting from the definition of t:
...
> t <- sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Set the column names to the gene names from one or the other of your conditions:
> colnames(t) <- rownames(immune)

Set the cutoff:
> cutoff <- 0.05

Filter for rows which meet the p-value cutoff criterium, and get their column names (gene names):
> colnames(t)[which(unlist(lapply(t["pval",], function(x) x < cutoff)))]
[1] "ABCC6"
...

Only one gene name is reported here, using your demo datasets. Adjust your cutoff and criterium logic, as needed. For example, I filter here for gene names where the p-value is less than the cutoff; however, from your question, you may be asking for the inverse, for whatever reason, so you could change < to >=.

I might suggest changing the variable name t to something else, as t() takes the transpose of the argument, if of matrix or data frame type. Using non-ambiguous variable names will make it easier for you to make changes or fix problems down the road.

From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:
pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]

Thank you, I added output of what you suggested in my post

Sorry, I have 2 matched matrices with matched samples and genes I want to find genes with p-value < 0.05 between 2 matrices but I could not do that

@FereshTeh forget your function f, there is no need, just run the commands in my answer. And check what you see in pvalueGenes object.

Thank you but your code itself giving error please have a look at error I added to my post

@FereshTeh The reason for error is statistics, read about it more here. Suggestion, when you get an error try to use the search, 99% of the time the error has been searched before and you will find the solution/explanation.

@zx8754, although names(x)[ x < 0.05 ] works, would it not be better to use names(x[ x < 0.05 ]), so each layer is a logical unit?

@RamRS Not sure what you mean by "logical unit". Yes both approaches work, but your version is about 1.14 times slower, tested on 100K vector.

Ah, I see. Never mind, then. By logical unit, I meant that x, x[x<0.05] and names(obj) are all self-contained, whereas names(x)[boolean_vector] is written in a way that the subsetting depends on understanding the statement as a whole.


text:
How do you query and explore ENCODE data?
I am looking for a modular way to query data from ENCODE.
For example, I would like to get CHiP-seq or similar tracks for a specific cell line. What's the proper way to do it?
Finally, is there an API to do it?

I don't know whether there is an API, but ENCODE's website does provide an interactive data matrix where you can filter data based on assay and sample type, place data sets in a "shopping cart", and then proceed to "checkout" to download the files of interest.

I'll add the fact that I've never used this website before today, so the fact that I was able to navigate it so easily is a credit to the group(s) that created and maintain the website!

You can just add /?format=json to any page to get the JSON output.
ENCODE REST API documentation: https://www.encodeproject.org/help/rest-api/
Example scripts: https://github.com/ENCODE-DCC/submission_sample_scripts

So let's say I want to get the bed file of this experiment using curl https://www.encodeproject.org/experiments/ENCSR000CKC/

If you look at the JSON for that experiment, there are multiple BED files listed there: https://www.encodeproject.org/experiments/ENCSR000CKC/?format=json

So let's say I want to get the bed file of this experiment using curl encodeproject.org/experiments/ENCSR000CKC
You can also search for files that belong to a given experiment, e.g. https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/

Then further select file properties (status, format, etc.) using the facets on the left.
This basically builds up a search URL for you, which is also possible to do manually. For example all BED files in ENCSR000CKC with GRCh38 assembly, returned as JSON: https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/&file_format=bed&assembly=GRCh38&format=json (See properties tab here https://www.encodeproject.org/profiles/file.json for all of the file fields you can use.)
You can also specify a field parameter to return only fields you are interested in. For example https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/&file_format=bed&assembly=GRCh38&field=s3_uri&field=cloud_metadata.url&field=dataset&field=file_format&format=json only returns the s3_uri, cloud_metadata.url, dataset, and file_format fields. These objects are returned in the @graph list:

Putting this all together to parse in Python:
>>> import requests
>>> import pandas as pd
>>> url = (
...     'https://www.encodeproject.org'
...     '/search/?type=File'
...     '&dataset=/experiments/ENCSR000CKC/'
...     '&file_format=bed'
...     '&assembly=GRCh38'
...     '&field=file_format'
...     '&field=assembly'
...     '&field=s3_uri'
...     '&field=cloud_metadata.url'
...     '&format=json'
... )
>>> r = requests.get(url)
>>> files = r.json()['@graph']
>>> df = pd.DataFrame(files)
>>> df['cloud_metadata'] = df['cloud_metadata'].apply(lambda x: x.get('url'))
>>> df[['@id', 'assembly', 'file_format', 's3_uri', 'cloud_metadata']]
                   @id assembly file_format                                                                                s3_uri                                                                                cloud_metadata
0  /files/ENCFF368QQU/   GRCh38         bed  s3://encode-public/2016/08/04/10493464-917f-4290-bb13-c5fec8319887/ENCFF368QQU.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/10493464-917f-4290-bb13-c5fec8319887/ENCFF368QQU.bed.gz
1  /files/ENCFF554NVI/   GRCh38         bed  s3://encode-public/2016/08/04/865f220e-7a47-427a-a7c1-9bdc049d4658/ENCFF554NVI.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/865f220e-7a47-427a-a7c1-9bdc049d4658/ENCFF554NVI.bed.gz
2  /files/ENCFF217CRP/   GRCh38         bed  s3://encode-public/2016/08/04/43aa693e-1581-4886-8ace-7e1609f45790/ENCFF217CRP.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/43aa693e-1581-4886-8ace-7e1609f45790/ENCFF217CRP.bed.gz

The s3_uri and cloud_metadata.url both point to the same thing. You can use the s3_uri to download using awscli or boto3. You can pass the https link to cURL.

There are two places:

encodeproject.org
screen.encodeproject.org

Details can be found in the ENCODE3 papers published in July 2020


text:
How do I re-name the headers of my Fasta file?
I appologize for asking this, but I really am really bad with regex...
Can someone help me transform the headers of my fasta files from this:
>Agaricus_chiangmaiensis|JF514531|SH174817.07FU|reps|k__Fungi;p__Basidiomycota;c__Agaricomycetes;o__Agaricales;f__Agaricaceae;g__Agaricus;s__Agaricus_chiangmaiensis

to this:
>JF514531|SH174817.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

seems like I could use the command
sed 's/__/:/g;s/;/,/g;s/.*|/|/g'

to substitute : for __, and , for ;.  I'm not entirely sure how to proceed from there though...
Thank you

Could you [edit] your question and explain the transformation you need please? You don't only want to replace __ with : and ; with ,. You also seem to be removing some text. Do you always want to remove whatever is before the first |? And what happened to the reps? And why is k__Fungi becoming tax=d:Fungi instead of k:Fungi? We need to know exactly what you are trying to do in order to give an answer that works for you.

Looks like you figured it out without more clarification - I'll try to be more concise when asking in the future -

If awk is good as well:
$ awk -v FS="|" -v OFS="|" '/>/{
    gsub(";",",", $5);
    gsub("__",":",$5);
    print ">"$2,$3";tax="$5"_"$3";";
    next;
}1' file.fa
>JF514531|SH174817.07FU;tax=k:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

This relies on the order of the parts separated by |.
I don't understand why k__Fungi becomes d:Fungi, so I ignored it for now.
Some explanation for what awkis doing.

| is defined as a field (column) separator for the input and output
If the line contains a > the code between { and } is executed.
the first gsub replaces ; with , in the 5th field
the second gsub replaces __ with : in the 5th field
then we print the new header line. The $2, $3 and $5 get replaced by the values of the second, third and 5th field of the input. , get replaced by the defined output field separator
with next we tell awk to go to the next input line and repeat the code
The 1 is in equivalent to print $0 for printing the whole input line. In our code this is only executed if the line is no header line.

I don't know awk that well. However, I find this is a really nice example and it makes me want to learn more awk!

Note that this won't print the sequences, only the header.

I just noticed this - how would you change this?

I edited my answer.

It seems like you want to do four things:

Remove all text up to and including the first |.
Now that that has been removed, take the string between what is now (after step 1) the first | and what is now the second |, and append it to the end of the header with a _.
Replace all __ with : and all ; with ,.
Change the (fixed?) string reps|k with tax=d.

If so, and assuming that you only have one occurrence of reps| in each header, you can do:
$ sed 's/>[^|]*|/>/;
       s/reps|k/tax=d/;
       s/__/:/g;
       s/;/,/g;
       s/|\([^|]*\)\(|.*\)/|\1\2_\1;/; ' file.fa
>JF514531|SH174817.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

This will:

s/>[^|]*|/>/ : replace the longest stretch of non-| characters after a > and until the first | with > alone.
s/reps|k/tax=d/; : replace the first occurrence of the string reps|k with tax=.
s/__/:/g; s/;/,/g : replace all (that's what the /g means) occurrences of __ with : and of ; with ,.
s/|\([^|]*\)\(|.*\)/|\1\2_\1;/ : the \(foo\) construct is called a capturing group and will capture foo and save it as \1, which can be used on the right hand side of the replacement operator. The next captured group will be \2 then \3 etc. Here, with |\([^|]*\) we are capturing the longest string of non| found after the 1st | (this happens after the first operation, so the 1st | is the one that was 3rd in the original) and then everything after the next | (|.*) until the end of the line. We then replace everything matched with a |, then whatever was matched by the first capture group (\1), then the second (\2), then we add a _ and repeat the \1.
Finally, if you want to make the changes to the original input file, use -i:
sed -i.bak 's/>[^|]*|/>/; s/reps|k/tax=d/; s/__/:/g; s/;/,/g; s/|\([^|]*\)\(|.*\)/|\1\2_\1;/; ' file.fa

That will save the original file as file.fa.bak. use -i alone without any argument if you don't need the backup.

This works really well, except for two things:  the command changes reps to tax for every other header, and there is still a | in front of tax instead of ;.


text:
Parallel in PLINK for linear association for SNP effects
I am doing linear association calculations for SNP effects, I see it take a lot of time to do, does anyone have any experience in parallelizing this job. I have used --threads <num of cores> but it says it only uses 1 thread
here is my code for the calculation in plink:
plink --bfile ukbb_dis --pheno phenotype_dis.pheno --mpheno 4 --allow-no-sex --covar phenotype_dis.pheno --covar-number 5-21 --missing-phenotype -10000000 --linear --adjust --out assoc_SNP_height_linear --threads 4

Here is the current output:
257682 MB RAM detected; reserving 128841 MB for main workspace.
1133273 variants loaded from .bim file.
284516 people (0 males, 0 females, 284516 ambiguous) loaded from .fam.
Ambiguous sex IDs written to assoc_SNP_height_linear.nosex .
284516 phenotype values present after --pheno.
Using 1 thread.
Warning: This run includes BLAS/LAPACK linear algebra operations which
currently disregard the --threads limit.  If this is problematic, you may want
to recompile against single-threaded BLAS/LAPACK.
--covar: 17 out of 21 covariates loaded.
Before main variant filters, 284516 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992737.
1133273 variants and 284516 people pass filters and QC.
Phenotype data is quantitative.
Writing linear model association results to
assoc_SNP_height_linear.assoc.linear ... 0%

Your help is really appreciated!

Try the latest plink 2.0 alpha build for this.  It’s an incomplete program, but its implementation of —linear is far better than v1.9’s.

Thank you! I have a little issue with the covariate-related. I have asked on Biostars here: https://www.biostars.org/p/359523/. Could you please have a look? Thank you!


text:
How do I prevent the FeatureHeatmap function from the Seurat package, from sorting my data groups in alphabetical order when plotting data?
How can I prevent a function from sorting my data groups (factors) in alphabetical order without affecting the integrity of the data? I am analysing single cell RNA sequencing data using Seurat 2.3.4. I want to be able to visualise the expression of two genes, Gene1 and Gene2 across cell clusters using the FeatureHeatmap() function from Seurat. I have assigned to the metadata slot of my Seurat object the factors or conditions being tested in the experiment. However on the resulting FeatureHeatmap plot, the data are sorted and grouped in alphabetical order to my disfavour. I have attempted a solution given on stackoverflow such as the one below:
conditions <- c(rep("WTNaClctrl", 1229),
               rep("WTBlu", 1328),
               rep("ILKOBlu", 349))

conditions <- as.data.frame(conditions)

seurat_object@meta.data$conditions <- factor(conditions$conditions, levels=unique(conditions$conditions))

When I use the FeatureHeatmap() function, the plot still sorts the data groups in alphabetical order:
FeatureHeatmap(seurat_object, features.plot = c("Gene1", "Gene2"),
               pt.size = 0.25,
               group.by = "conditions",
               max.exp = 3)

This plot gets produced:
I would like groupings not to appear in alphabetical order as such:

As previously suggested StackExchange, I have tried to upgrade my Seurat object to Seurat V3:
UpdateSeuratObject(Seurat_object):
Object representation is consistent with the most current Seurat version.

devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0') : package ‘curl’ successfully unpacked and MD5 sums checked:

Error: (converted from warning) cannot remove prior installation of package ‘curl’.

Thank you in advance for any suggestions.

The error you got is from not having the rights to update curl package. Maybe it was installed by the administrator of the computer (sudo?). The plotting side of seurat seems to be done with ggplot2, which orders the categoric values based on the levels of the factors. Reorder the levels and it should go as you wish (I'd recommend forcats)

@IIrs Please I am not clear about the point in the code where I should apply forcats.

forcats is a package designed to make it easier to work with factors, it provides some functions to reorder the levels of the factors. I don't know seurat good enough to suggest you where you should use it.

Possible duplicate of How to set the position of groups in a Seurat object on a FeatureHeatmap plot

@TimStuart: I think in the other question, the alphabetical impact is not highlighted. So I thought of making it a new question.

As I said in your previous question, you can't change the order using Seurat v2. The groups are coerced to a character vector which is then ordered alphabetically by ggplot2. You can either use Seurat 3 and reorder the factors, or extract the data from the object and plot it yourself

@TimStuart: In the question I also highlighted my difficulty in automatically downloading the Seurat V3 package due to admin rights issues. Could I please have a link where I could download the package and load it manually?

See the github page here: https://github.com/satijalab/seurat/tree/release/3.0 and the website here: https://satijalab.org/seurat/install.html#prerelease

I had this same issue. I reordered the levels and that worked for most functions, but would switch back with the FeatureHeatmap. Turns out there is a line in the function that changes it to a character.
I ended up making a "new" function and just changed that one line from:
data.plot$ident <- as.character(x = object@ident)

to:
data.plot$ident <- object@ident

Once that was fixed, the order that I assigned worked out just fine.


text:
Calculate molecular mass from chemical formula C9H15N4O8P with Python?
Can I calculate the molecular mass from a formula like this C9H15N4O8P with Python e.g. using rdkit? I think rdkit.Chem.Descriptors.MolWt needs a molecule object. Though, to generate that I would need e.g. a SMILES string. Is there any other python module that can do the trick?
I found this: http://verahill.blogspot.com/2012/10/molecular-weight-calculator-in-python.html

I tried to adapt this approach within rdKit to calculate exact mass (for mass spectrometry), but was not successful.
This library (https://pypi.org/project/molmass/) appears to be able to do this:
from molmass import Formula
f = Formula('C7H11N3O2')
f.formula
f.isotope.mass
'169.08512660696'

As you mentioned, rdKit does not accept a chemical formula as input. The main reason for this are probably isomers. You could however use rdKit to initiate all the atoms, get their mass and add them together.
import re
from rdkit import Chem

def get_mass(formula):

    parts = re.findall("[A-Z][a-z]?|[0-9]+", formula)
    mass = 0

    for index in range(len(parts)):
        if parts[index].isnumeric():
            continue

        atom = Chem.Atom(parts[index])
        multiplier = int(parts[index + 1]) if len(parts) > index + 1 and parts[index + 1].isnumeric() else 1
        mass += atom.GetMass() * multiplier

    return mass

Then...
>>> get_mass("C9H15N4O8P")
338.21299999999997


text:
structural variants: why doesn't AS + AP = SU?
I have run smoove then duphold on 46 individuals, then pasted the duphold output back into a squared vcf with all 46 samples.
I have been using SU/AC as one of my filters, looking for variants that have a minimum threshold of support per alternate allele called. Now I'm confused by some of the variants that are getting filtered out, and wonder why SU (number of pieces of evidence supporting the variant) is not calculated as I had expected.
For one variant, if I sum the AS values (split-read observation counts) and the AP values (paired-end observation counts) over all 46 individuals, I get AS=24 and AP=15, so 39 total reads supporting the variant in the population. However, SU is only 14 (PE=5 + SR=9).
Can anyone tell me why the values for SU, PE, SR are so much lower than the total AS & AP counts? Thank you!

We welcome this question here, but you may get a quicker and more thorough response posting to the smoove issue tracker.

Thanks - it's now posted here: https://github.com/brentp/smoove/issues/53


text:
Recommendation for a binary dataset for Computational Biology open access (as in a dataset with 1 and 0 that can be used to apply ML techniques?
I have a couple of modified ML methods that I want to test with biological binary data - specifically binary gene expression data - e.g. https://eprints.soton.ac.uk/69359/1/Thesis.pdf that might include a variable that could be used to be the target variable in ML classification. What I mean by binary data is the features are all either 0 or 1, so I expect a large dataset e.g this thesis. I am looking for any suggestions.

Please [edit] your question and make it a bit more specific What kind of "biological data" do you need? How large a dataset? What do you mean by "binary"? Are you referring to the format of the data or do you want something that can only have two states?

Kaggle has a cats v dogs dataset if you want to try your hand at biological image classification: https://www.kaggle.com/c/dogs-vs-cats
The DREAM challenge competitions are usually specifically oriented to solving biological classification or learning problems with ML: http://dreamchallenges.org/


text:
RNA-seq gene/transcript read counts database for Mouse
Is there any RNAseq gene/transcript read count database for the mouse? I already know about ARCHS4, looking for some other source.
Thank you

You can sometimes get counts from GEO, though that's probably not what you're looking for.

EBI has an expression atlas, including mouse: https://www.ebi.ac.uk/gxa/experiments?species=mus+musculus
If you click on one of the individual experiments you can download tables of TPM or RPKM per gene.
There are currently about 1000 individual experiments in different tissues etc.


text:
Python - Finding a motif - input: a txt file with 10 sequences and 10 motifs
When I run my BruteForce function with only one input it works and the result is correct.
def BruteForce(s, t):
    occurrences = []
    for i in range(len(s)-len(t)+1): # loop over alignment
        match = True
        for j in range(len(t)): # loop over characters
            if s[i+j] != t[j]:  # compare characters
                match = False   # mismatch
                break
        if match:   # allchars matched
            occurrences.append(i)

    print(occurrences)

t = 'CAACTTCCA'

s = 'GACAACTTCCAACTTCCAACTTCCCGTCCCAACTTCACAACTTCGGCCCAACTTCCATGCAACTTCACCATCAACTTCGCTCGAAGCTGCCTTCCACTCCAACTTCACAACTTCCTCAACTTCCTCACCAACTTCAGCAACTTCTCTAGGGCCAACTTCCAACTTCTCAACTTCTCAACTTCCAACTTCCGACAACTTCTCCTGGCAACTTCCAACTTCCAACTTCAATACAACTTCGCAGACAACTTCCGCAACTTCGAACAACTTCCAACTTCCCCAACTTCCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCCCAACTTCAGATAGCAACTTCGATCTTACACAACTTCACGCAACTTCTCCAACTTCCAACTTCTGTGCAACTTCTCTGAACAACTTCCTCAACTTCCAACTTCGCAACTTCCCCAACTTCCTCAACTTCATGCAACTTCGAGGCAACTTCCCAACTTCGCAACTTCCTATTCCCAACTTCTGTGGCAACTTCTCAACTTCTGGACAACTTCTATGCCCAACTTCACAACTTCCCCAACTTCTTTACAACTTCGACAACTTCATCAACTTCTAGTCAACTTCTGGTCCAACTTCCAACTTCCCCAACTTCCAAAGTGCCGCAACTTCGTAACAACTTCACGCGCTCAACTTCAACCAACTTCTTTTCCCGCAACTTCGCAACTTCACAACTTCTAATCAACTTCCAACTTCGGATCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCTCCAGGGACAACTTCAAGTACAACTTCCAACTTCGCAACTTCACAACTTCCCAACTTCGCAACTTCTACACGCAACTTCCAACTTCTGGTCCCAACTTCATCAACTTCAGTCAACTTC'

BruteForce(s, t)

[2, 9, 48, 152, 175, 205, 212, 261, 277, 284, 300, 307, 370, 410, 607, 623, 717, 735, 751, 758, 792, 844]

But when I use it in the following loop it returns nothing! I have used the print() to see if it reads the strings and motifs from the input_2.txt file and it does but it returns an empty list of occurrences [].
with open('input_2.txt', 'r') as input_2:
    nt = input_2.readline() #nt is the no of strings+motifs in 1st line
    nt = int(nt) #in this case it is 10

    for i in range(nt):
        s = input_2.readline()
        print(s)
        t = input_2.readline()
        print(t)
        BruteForce(s, t)

Could you please let me know what changes are needed here?
Many thanks,
Beeta
input_2.txt:
10

GACAACTTCCAACTTCCAACTTCCCGTCCCAACTTCACAACTTCGGCCCAACTTCCATGCAACTTCACCATCAACTTCGCTCGAAGCTGCCTTCCACTCCAACTTCACAACTTCCTCAACTTCCTCACCAACTTCAGCAACTTCTCTAGGGCCAACTTCCAACTTCTCAACTTCTCAACTTCCAACTTCCGACAACTTCTCCTGGCAACTTCCAACTTCCAACTTCAATACAACTTCGCAGACAACTTCCGCAACTTCGAACAACTTCCAACTTCCCCAACTTCCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCCCAACTTCAGATAGCAACTTCGATCTTACACAACTTCACGCAACTTCTCCAACTTCCAACTTCTGTGCAACTTCTCTGAACAACTTCCTCAACTTCCAACTTCGCAACTTCCCCAACTTCCTCAACTTCATGCAACTTCGAGGCAACTTCCCAACTTCGCAACTTCCTATTCCCAACTTCTGTGGCAACTTCTCAACTTCTGGACAACTTCTATGCCCAACTTCACAACTTCCCCAACTTCTTTACAACTTCGACAACTTCATCAACTTCTAGTCAACTTCTGGTCCAACTTCCAACTTCCCCAACTTCCAAAGTGCCGCAACTTCGTAACAACTTCACGCGCTCAACTTCAACCAACTTCTTTTCCCGCAACTTCGCAACTTCACAACTTCTAATCAACTTCCAACTTCGGATCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCTCCAGGGACAACTTCAAGTACAACTTCCAACTTCGCAACTTCACAACTTCCCAACTTCGCAACTTCTACACGCAACTTCCAACTTCTGGTCCCAACTTCATCAACTTCAGTCAACTTC

CAACTTCCA

ATTGTGTATCGCTATGTATCGTGTATCGGATTTTGTATCGTGTATCGGTGTATCGTGTATCGTGTATCGTATCACTGTATCGTTGTATCGTAGCGTTGTATCGTAATGTATCGCTCTGTATCGTGTATCGGGTTTGTATCGATGTGTATCGCTGTATCGGTGTATCGGTGTATCGCTTGTATCGACAGCGCTTGTATCGTGTATCGACCTGTATCGGTGTATCGTGTATCGAATGTATCGTTGTATCGAATTGTGTATCGTGTATCGTGTATCGTGTATCGTATGTATCGACTGTATCGCTGTATCGTAGCCTGTATCGTGTATCGGTGTATCGGTGTATCGGCGTGTATCGAATGTATCGTTGTATCGCTGTATCGTGTGCTGTGTATCGGTGTATCGACCCGTTGTATCGTGTATCGTGGGTAAGTGTATCGTTGTATCGTAGTGTATCGTGTATCGTGTATCGTCATGGTATTGTATCGTTATAGCTGTATCGCTCGGCTGTATCGTATGTATCGTGTATCGTGTATCGCATGTATCGGTGTATCGCTGTATCGACCATGTATCGTGTATCGGTGTATCGATGTATCGCTTCCATTAGAAATGTATCGTATGTATCGTGTATCGCTGTATCGTTTGTATCGCATGTATCGATTGTGTATCGTGTATCGTTGTATCGTGTATCGTGTATCGTTGTATCGTAATAACGATGTATCGAATGATGTATCGTTGTGATGTATCGTGTATCGACTGTATCGATGTATCGGGTTGTATCGTGTATCGTGTATCGTAGATGTATCGAGAGCCATTGTATCGCCTTGTATCGCGTGTATCGGTGTGTATCGTTTGTATCGTGTTTGAATGTATCGCTGTATCG

TGTATCGTG

TTTCAGCGTGTTTGGCTTTCAGCACTACTTTCAGCGGTTTCAGCTTTTTCAGCTTTTCAGCTATTTTCAGCTTTCAGCCATCTTTCAGCCCAGGTTTCAGCTGTTTCAGCTTTCAGCGAGTTTCAGCTTTCAGCGCGGGGATTTCAGCGTATTTCAGCGACATTTCAGCAGTTTCAGCTGGAGTGAAGCCGTTTCAGCGCTTTTCAGCACCATTTCAGCTAGTAGGTTTCAGCTTTTCAGCCTTTTCAGCTTTCAGCTATTTCAGCTCTAGAAGTCGTTTCAGCTTTCAGCTTCCTTTTCAGCTTTCAGCTTTCAGCCTTTCAGCTTTCAGCAAACTTTCAGCGATTTCAGCTTTCAGCTTTCAGCTTTCAGCTTTCAGCATTTCAGCTTTCAGCCGTTTTCAGCGACTCTTTCAGCGTTTCAGCTTTCAGCTGATCGTTTTCAGCTTTTCAGCCGGTTTCAGCGAAAGTTGGTCTTTTCAGCCAAATTTTCAGCTTTCAGCTTTCAGCGTTTCAGCTTTTCAGCCATTTCAGCCTATTTTCAGCAATCTTTCAGCAATTTTCAGCGCAATTTCAGCCAAAATTTCAGCATTTCAGCTTTCAGCGTAGTTTCAGCTTTCAGCTGAGCCTTTCAGCTTTTCAGCTCTCTTTTCAGCTTTCAGCGATTTCAGCCTTTTCAGCGCTTTCAGCACCACGCCTTTCAGCTTTCAGCTTTTCAGCCTTTCAGCCTATTTTCAGCGGTTTCAGCGTTTCAGCTTTCAGCACTTTCAGCTTTCAGCTTTTCAGCGGAAGGTTTTCAGCA

TTTCAGCTT

AGCATCCGACTGCATCCGTGCATCCGACGGCATCCGGCATCCGCAGCATCCGCACAGCCGCATCCGAGCCCGCATCCGAAGCATCCGGCATCCGGGCATCCGGCCGCATCCGGCATCCGGGCATCCGTGCATCCGGCATCCGTGTGCATCCGGCATCCGGGGTGCTTGCATCCGCGCATCCGTGCATCCGCGCATCCGGCATCCGCAGTCCGCATCCGCGCATCCGGCATCCGAGCATCCGTATTCGCATCCGTCGAGCATCCGTACGTCGCATCCGTGCATCCGGGAGAGGCATCCGGGCATCCGCGCATCCGGTGGACTATAACGCGCATCCGGCGAGCATCCGCCGCGGCATCCGTGCATCCGGTACGGCATCCGGGCTGAGGCATCCGCTCCCGGGGCATCCGCGGCGCATCCGCACGCGCATCCGGCATCCGAGCATCCGTTGCATCCGGAGCATCCGTAGCATCCGATAGCATCCGCAAGCTTGCATCCGGCATCCGCTACAGCATCCGGTGGCATCCGAAGCATCCGCGCATCCGCGGATTGCATCCGCACACGGCATCCGTAGCATCCGTGCATCCGGAGTGCATCCGCCGGGCGCATCCGCTGCATCCGCAGCGCATCCGAGCATCCGATCTTGCATCCGGCATCCGGGAATGCATCCGGGACTGCATCCGGTCTTAAGGGTGCATCCGAATGCATCCGCTGTAGCATCCGGCATCCGAGTAAGCATCCGAGTTCTGCATCCGGAAGCAGCATCCGGCATCCGGACACCAGCATCCGCGCATCCGGGGCGAGCATCCGAGCATCCGGCATCCGGGGCAAGTGGCATCCGGCATCCGTCGTGCATCCGGGCATCCGAA

GCATCCGGC

GACGAGCTGACGAGCGTGACGAGCAGGCGACGAGCGACGAGCATGCGGTGACGAGCGGACGAGCGACGACGAGCGCGGACGAGCGACGAGCGGACGAGCTCCGACGAGCGACGAGCACGACGAGCGCCCGACGAGCGACGAGCATGATGACGAGCGACGAGCTTAGACGAGCATGCTGACGAGCTGACGAGCTGACGAGCTCGTACATCGACGAGCGGACGAGCAGCCCGATAAGCCTTCGACGAGCTGACGAGCCACGACGAGCGACGAGCGACGAGCGGACGAGCGGGACGAGCGTCGGACGAGCGGACGAGCAGGACGAGCACCTCAATCGACGAGCTGACGAGCGGACGAGCCGACGAGCTGAACTGAAGGACGAGCCTACGTTCTAACGTGCCGTCACTGACGAGCGACGAGCGGACGAGCGATGGACGAGCTGACGAGCGCAGGACGAGCGAGAAGGCTGACGAGCAGACGAGCTGACGAGCTCGACGAGCAGACGAGCGGACGAGCTTTGATAAGACGAGCACGTCGACGAGCCTCGTGACGAGCTGACGAGCGACGAGCTAGACGAGCGACGAGCTGACGAGCGACGAGCACTAAACGCGACGAGCTCGACGACGAGCATGGACGAGCGACGAGCTAGGACGAGCGGCGACGAGCAGACGAGCGACGAGCTGACGAGCTTATAGACGAGCCTGACGACGAGCCAAGACGAGCGACGAGCGGACGAGCGACGAGCACGACGAGCCGACGAGCCGCGGACGAGCCGTAGACGAGCCAATCATTAAGACGAGCCGACGAGCCATTTGGGGACGAGCCTCCTCGACGAGCCATAAGACGAGCTGACGAGCCATGACGAGCATGCCCGACGAGC

GACGAGCGA

GAGGACCCCGGACCCCTGGACCCCAGGACCCCGACGGTGGGGCGGACCCCGAGAAGGACCCCTGGACCCCCGCAGGACCCCTTTATCGGACCCCGGACCCCGGACCCCCGGACCCCGGCTGGACCCCGTCGTAAGGACCCCGAATCGGACCCCTAGGACCCCAGGGACCCCTCCCCGGTTGGACCCCCTGGACCCCCTTGAGAGGGACCCCGGACCCCACGCCGTGCTTAAGGACCCCACTATGGACCCCATACGGACCCCGGACCCCGATCAGAGCGACCAGGACCCCCGGCCTGGACCCCTCGGACCCCAGAAGGACCCCTCGGACCCCTGGACCCCGGACCCCACAGGGACCCCGGACCCCTAGCCGCGGTGGACCCCCGGACCCCCGCAGATGGGGACCCCTCATGCGGACCCCCTGACGGACCCCTGGACCCCCGGGACCCCCAGGGACCCCATTTCGAGGACCCCATGGGGACCCCAAGCTGGACCCCGGACCCCTGGGACCCCCGGACCCCGGACCCCGGACCCCATAGGACCCCCGTTTGTTGCCATGGACCCCTTGGGACCCCGAGTCGGACCCCGGACCCCCAGGACCCCACAGGACCCCGGGACCCCGGACCCCATCGATCGGACCCCGGACCCCGGACCCCTACGGACCCCGCTAAGGGACCCCGGACCCCGTGGACCCCCACGGACCCCTAGGACCCCGGGACCCCGGACCCCGGACCCCAGTTGGACCCCCCTTAGGACCCCAGGACCCCACATGAGGACCCCGGGACCCCGGGACCCCTGGACCCCTGGACCCC

GGACCCCGG

AAACCTTCGGACCTTCGTACCTTCGACCTTCGTGAAACCTTCGGACCTTCGTGACCTTCGGGTGACCTTCGGGTACTCACCTTCGGCACCTTCGTATACCTTCGAGACCTTCGCTGGACTGTAAACCTTCGATCGTTATTTTCTACCTTCGCACACCTTCGACCTTCGAGCCGAACCTTCGGACTGGCCTCACCTTCGATCACCTTCGACCTTCGAACCTTCGACCTTCGACCATACCTTCGACCTTCGGACCTTCGGACCTTCGCACCTTCGAACCTTCGCCACCTTCGACCTTCGAACCTTCGGTGGGCCACCTTCGTTACCTTCGAACCTTCGGGCACCTTCGACCTTCGGCTTACACCTTCGGACCTTCGATGACTGACCTTCGAAACCTTCGACCTTCGACCTTCGATTCCACACCTTCGACCTTCGACACCTTCGACCTTCGCGTTCACCTTCGACCTTCGACCTTCGGGCCAATACCTTCGACCTTCGGCCTAATACCAACCTTCGAACCTTCGAACCTTCGGGTGACCTTCGGACCTTCGCGCACACCTTCGCACTTACCTTCGACCTTCGGTGACCTTCGACCTTCGTGCACCTTCGCGTCTCACCTTCGCACCTTCGTTACCTTCGAACCTTCGACCTTCGACAACCTTCGTACCTTCGCTACCTTCGGGAACCTTCGATGTCCTACCTTCGACCTTCGCGACCTTCGACCTTCGTAACCTTCGAACCTTCGACCACCTTCGAACCTTCGAAGGAAGGACCTTCGAATCATTACCTTCGTTACCTTCGTAGTACCTTCGCCGACCTTCGGACCTTAACCTTCGTGA

ACCTTCGAC

AAGTGGTATTCCCAAAAAGATATTCCCGTTATTCCCGCTATAAATATTCCCAACTGTTTATTCCCACTTATTCCCGTATTCCCGTATTCCCCGGGCCTTTTATTCCCCATTCTATTCCCGTATTCCCATATTCCCACTTATTCCCCTATTCCCATATTCCCTATTCCCTATTCCCGTAAGTATTCCCCCCCTTATTCCCCCATTTGTCCATATTCCCTCAATATTCCCTATTCCCCTATTCCCGCTTCCCCAGTCTATCCCAAATTGTATTCCCAATATTCCCCTTTATTCCCTTTATTCCCGTATTCCCTATTCCCGTGGCTATTCCCCCGACGTTATTCCCTCTATATTCCCGTCTCAAAGTATTCCCTTATTCCCAGTATTCCCATTCGCCTGGTATTCCCGTATTCCCGGTATTCCCGTTTATTCCCATCTATTCCCATTATTCCCTGGTATTCCCGCTATTCCCGCAACGGTATTCCCTATTCCCTATTCCCGTATTCCCAATATTCCCTTCCTATTCCCGTATTCCCTAGCTATTCCCATATTCCCGATATTCCCTATTCCCGTATTCCCTATTCCCTTATTCCCTATTCCCTAGTTATTCCCAAACAGGCTATTCCCTATTCCCTATTCCCATATTCCCTTATTCCCTACGTATTCCCTCTTAAAAATATTATTCCCCTATTCCCTATTCCCGTATTCCCCGAAAGTATTCCCCCATATTCCCAATTATATTCCCGGTATTCCCCCGTATTCCCTATTCCCTTATTCCCATATTCCCATATTCCCAAAGTTTATTCCCGTATTCCCCACTATTCCCATAGATTATTCCCTATTCCCTATTCCCTTATTCCCCTATTCCCTATTCCCTTATTCCCGTATTCCCAGTTTATTCCC

TATTCCCTA

GAAAGGATAAAGGATGTCAAAGGATGCAGATATAAAGGATAAAGGATACAAAGGATTAAGTATGTCCGAAAAAGGATAAAGGATTAAAAGGATCGCTGAACCTTACACAAAGGATAGTGAACAAAGGATTCAATAAAAGGATAAAGGATCAAAGGATCAAAAGGATACAAAGGATGCCGATGAAAGGATAAAGGATCTAAAGGATAAAGGATGAGAAAGGATTGGAAAGGATTATGAGATCAAAGGATAAAGGATGGTGTAAAAGCTAAAGGATTCGGCAAAGGATAAAGGATTAAAGGATAAAAAGGATAAAGGATGAAAGGATCCGCAGGGACCAGCAAAAGGATAAAGGATCGAATGGGTAAGAAAGGATCAAAGGATGAAAAAGGATTACTAAAGGATAAAGGATCCTGAAAGGATTACAAAGGATCTTAAAAGGATTCGGAAAGGATCCATAGGAAAAGGATAAAGGATCGCGAAAGGATAAAGGATTAAAGGATAAAGGATAAAAGGATTCAAAGGATAAAGGATAGACAAGGGAGAAAGGATGCAAAGGATGAAAGGATAAAGGATAAAGGATAAAGGATTAGGTTAAAGGATCGAAAGGATCCAAAGAGAGCGAAAAGGATGAGGACGAAAGGATCAAAGGATCCAAAGGATCAAAGGATAAAGGATTGCGATGGAAAGGATGTCAAAGGATCACCAAAGGATAAAGGATAATAAAGGATAAAGGATAAAGGATCAAAGGATTTGAAAAGGATAAAGGATCACGAAAGGATAAAAGGATTGGCAAAGGATGAAAAGGATGAAAGGATAAGCCCTCCAAAGGAT

AAAGGATAA

ATTCAAATACTTCAAATAGTCAAATATCAAATATGCTTCAAATATCAAATATCGGTCAAATAAGTCAAATAAGTCAAATATCAAATATCAAATACTCAAATATCAAATAGTGTCAAATACGAATTGGGTCAAATATCAAATATGTGTTCAAATATTCTTCAAATACTGGACACTCAAATAGGAGTCAAATATCAAATATCAAATACGTGAAGTGCTTGTCAAATATTTCAAATACTTTCAAATATTCAAATACTCAAATATGTTCAAATATCAAATATCAAATATTTCAAATATCAAATATCAAATAGTCCTCAAATAGCAAACCAGTCAAATATCAAATATGTCAAATATGCTCACGGCAACCTCAAATATCAAATATCAAATAGATGTCAAATAAGTTTCAAATATCAAATATCAAATATCAAATATCAAATATTCAAATATCAAATAGCGGGCTCAAATAAGCTCAAATACGTCAAATAGGGGGTCAAATAGATCAAATAGTCAAATATTCAAATACATCAAATAGTCAAATACAAGAACCACCGAGATCTCAAATAATCAAATATCAAATATGGTCAAATAAATATTCAAATAGGTCAAATAAGATCAAATATCAAATAAGTCGTCCATCAAATAGTCAAATAGCTCAAATAACCTCAAATATCAAATATTCAAATACCGGTCATCAAATAGGACAAATCAAATAGTCAAATAAGATCCTCTCAAATAATTCAAATAGCTGTTCAAATACTCAAATATCAAATAGTCAAATATCAAATATTCAAATATCAAATACTTCAAATATGTTCAAATAATCAAATACTCAAATATCAAATAATCAAATAATACTTCAAATATACCAAACGCTCAAATATTAGTTGGATCAAATATCTTCAAATATCAAATAA

TCAAATATC

I tried to fix the formatting, but please check and make sure I didn't screw anything up. For help with formatting, please see formatting tools.

why not to use kmp? https://it.wikipedia.org/wiki/Algoritmo_di_Knuth-Morris-Pratt

Let's start by doing some explaining about why your code doesn't work.
with open('input_2.txt', 'r') as input_2:
    nt = input_2.readline() #This is reading the first line of the file (which is 10)
    nt = int(nt) # this is storing "10" as an int

    for i in range(nt): # This is the same as "for i in range(10):" which just returns nums 0-9, does not access lines of your file
        s = input_2.readline() # this is reading the second line which is blank
        print(s) # this will print a blank line 10 blank lines
        t = input_2.readline() # this is reading the rest of the file
        print(t)
        BruteForce(s, t) # this is calling your function on 10 blank lines and the restof your file, not each seq/motif

I made some modifications that will help you achieve your goal I think:
def BruteForce(s, t):
    occurrences = []
    for i in range(len(s)-len(t)+1): # loop over alignment
        match = True
        for j in range(len(t)): # loop over characters
            if s[i+j] != t[j]:  # compare characters
                match = False   # mismatch
                break
        if match:   # allchars matched
            occurrences.append(i)

    print(occurrences)

def group(lst, n):
    """ This function will be used to group your sequences and motifs later """
    for i in range(0, len(lst), n):
        val = lst[i:i+n]
        if len(val) == n:
            yield tuple(val)

with open('input_2.txt') as input_2: #standard open in read mode,
    lines = input_2.readlines()[2:] # skip 10 and the blank space bc we don't care about them
    split_line = [line.strip('\n') for line in lines if line!='\n'] #this will store lines that aren't blank into a list and strip new_line chars '\n'
    seq_tuple = list(group(split_line, 2)) # this will use group func above to group the list into tuples of 2's which are sequence and motif.

for x in seq_tuple: # this accesses the list of tuples (seq_tuple)
    BruteForce(x[0], x[1]) # x[0] in each tuple is seq, or BruteForce(s), x[1] in each tuple is motif, or BruteForce(t)

Running this will output Each motifs start position in each sequence:
[2, 9, 48, 152, 175, 205, 212, 261, 277, 284, 300, 307, 370, 410, 607, 623, 717, 735,
751, 758, 792, 844]
[14, 33, 48, 55, 116, 192, 217, 253, 260, 267, 312, 370, 405, 412, 445, 452, 511,
518, 561, 613, 657, 672, 679, 735, 769, 776, 847]
[37, 46, 64, 103, 120, 226, 243, 277, 284, 296, 303, 318, 345, 352, 359, 366, 381, 418, 439, 488, 495, 510, 590, 608, 628, 648, 697, 704, 748, 764, 771]
[29, 80, 95, 105, 128, 145, 191, 219, 328, 423, 489, 642, 713, 758, 808, 831]
[28, 57, 77, 102, 129, 148, 258, 265, 404, 419, 449, 553, 569, 584, 629, 664, 713, 728]
[2, 87, 94, 109, 204, 254, 332, 350, 487, 511, 518, 577, 603, 611]
[18, 154, 201, 216, 223, 235, 283, 339, 390, 397, 418, 425, 434, 453, 460, 481, 565, 582, 637, 644, 695, 711, 735]
[154, 161, 221, 303, 476, 483, 526, 554, 569, 584, 591, 617, 624, 647, 685, 754, 829, 836, 859]
[1, 33, 70, 135, 182, 198, 241, 279, 294, 303, 340, 395, 460, 478, 493, 500, 517, 559, 566, 573, 661, 704, 711, 721, 728, 754, 772, 808]
[19, 37, 44, 73, 80, 95, 128, 184, 191, 263, 270, 286, 293, 327, 364, 371, 400, 407, 414, 421, 436, 561, 609, 659, 754, 769, 784, 826, 886, 896]

Yes, yes, yes! You are right (and I feel very dumb), I should have paid attention to removing \n. All the times I spent on this problem.... I am thinking I could have changed the input file to something without any blank line and see how different the output would be. Thank you SO much for your time and your courtesy. I truly appreciate it. You have no idea how much this pipsqueak problem has been bugging me since yesterday. Thank you!

No problem! I’m happy to help, and welcome to the site! Going forward most bioinformatics file formats will utilize ‘\n’, ‘\t’, or ‘,’ as file separators in some form or another so it is good good to remember quick ways (like the list comprehensions above) to quickly manage those in your files. Thanks for your comment, it was very kind!

This looks hard work. In Perl the function "index" and "rindex" as well as "substr" will perform this without the incremental loop. There must be an equivalent in Python?
I've not tried this, but at a guess it would be something like,
my @array;
my $target = 'AGATAGAGAGATAGAGAGATAGAGAGATAGAGAGATAGAG';
my $query = 'AGAGA';
while(length($target) > 4){
    my $location = index ($query, $target);
    push @array, $location;
    $target = substr ($location, length($target), $target);
}
print join("\s", @array), "\n";

Thank you! There is the find() method which returns ONLY the first occurrence (not all of them). But I just found out that I haven't tried Regex in python.

Perl is the same. The idea is that if you use it recursively with a "string splicing" algorithm you can get all your motifs easily and low computuational overheads. Please check "find()" for something like rfind() .. which does the same thing from the oppositive side of the string. That can be useful because when find() = rfind() then all the motifs in the string have been identified.


text:
Why do I get so many insertions from Minimap2 on my Nanopore WGS?
I'm a starting my analysis on nanopore whole genome sequencing. I start my analysis from this popular Github.
The sample I downloaded was WGS for NA12878, so I would assume it's alignment to GRCh38 shouldn't be that bad.
But ... I'm getting lot's of insertions...

I'm not sure it was caused by the sequencing platform, the MiniMap aligner or my bugs.
I supply all my commands, starting from sample download.
wget https://s3.amazonaws.com/nanopore-human-wgs/rel5-guppy-0.3.0-chunk10k.fastq.gz
gunzip rel5-guppy-0.3.0-chunk10k.fastq.gz
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -d GRCh38_full_analysis_set_plus_decoy_hla.fa.index GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -ax map-ont -t 20 GRCh38_full_analysis_set_plus_decoy_hla.fa.index rel5-guppy-0.3.0-chunk10k.fastq | samtools view -bS - | samtools sort -@ 5 > rel5-guppy-0.3.0-chunk10k.bam
samtools index rel5-guppy-0.3.0-chunk10k.bam

Insertions and deletions are the dominant error mode of long read sequencing, including nanopore sequencing. What you see is not unexpected. Things may have improved by now if you would download the raw fast5 data and repeat the basecalling. There is no need to gunzip the fastq.gz prior to alignment. Your commands for alignment look alright, except that (if you are using a recent version of samtools) you don't really need samtools view and can pipe directly to samtools sort, using -o to specify the output format:
minimap2 -ax map-ont -t 20 <index> <fastq> | samtools sort -@5 -o alignment.bam

In a recent version of igv you can right-click the alignment and chose a "Quick consensus mode" (or something like that), which will make your alignments at least a bit prettier. Note how the consensus sequence (top of the graph) is quite clean - while individual reads may be noisy this is mostly randomly distributed and as such ironed out in consensus.


text:
Calculate average parental genotype
I am interested in creating an averaged value for genotype across two parents:
E.g. Input (three unrelated couples, 1, 2 and 3):
Individual    SNP1    SNP2    SNP3    SNP4 ...
Father_1      1       0       1       2    ...
Mother_1      0       0       2       1    ...
Father_2      1       0       1       1    ...
Mother_2      2       1       0       1    ...
Father_3      1       2       0       1    ...
Mother_3      1       .       0       2    ...

E.g. Output (the midparental genotypes for couple 1, 2 and 3):
Individual    SNP1    SNP2    SNP3    SNP4 ...
Midparent_1   0.5     0       1.5     1.5  ...
Midparent_2   1.5     0.5     0.5     1    ...
Midparent_3   1       .       0       1.5  ...

I can think of a way to do this using e.g. a PLINK .raw file and manipulating that in R, but am concerned that would become memory- and storage-heavy when assessing 100000s of SNPs and 1000s of pairs.
Does anyone know of an existing tool that can do this?
Thanks

Hi there, maybe it's just my ignorance, but what "average parental genotype" should represent? Two homozygotes (0 2) and two heterozygotes (1 1) have the same average (1) but produce quite different genotypes (0 1 2 in one case 1 in the other). Performance wise, would guess that a small python script could do this quite fast. 100,000 SNPs is not that much.

Welcome to the site Joni. Could you show us how would you do it using PLINK and R? This way we could improve the code or suggest faster/more efficient alternatives knowing the base line. To do so, remember you can [edit] your question.

Question is also posted on biostars.

Example pseudocode of how this could be achieved in PLINK and R. 1) Split PLINK binary file to couples and offspring. 2) Convert PLINK couple binary file to a .raw file with --recode A. 3) Sort by FID (assuming FID identifies couples). 4) Create new PLINK midparent.raw file in R by averaging each SNP across the two lines. 5) Create additional files needed to convert between PLINK .raw and PLINK binary. 6) Create new PLINK binary with midparent.raw. 7) Merge with offspring binary.

Beside the correct question by @Kamil, this sounds like a job for datamash.
First let's get the number of columns in the file
$ head -n1 genotypes.txt| awk -v FS="\t" '{print NF}'

Now lets calculate the average per column with datamash
$ datamash -H mean 2-<n> < genotypes.txt > mean.txt

Replace <n> with the number of columns we received with command above.
For your example the output will look like this:
mean(SNP1)  mean(SNP2)  mean(SNP3)  mean(SNP4)
0,5 0   1,5 1,5

One can remove the mean with sed
$ sed -i -e 's/mean(//g' -e 's/)//g' mean.txt

Sorry, this is my fault for not providing good enough example data - there are many pairs of rows in the file - will edit the example data to make this clear.

These are unusual parents because they all appear haploid. Anyway, if I was doing this I would calculate the allele frequency per progeny such as makefreq in adegenet in R-cran (just my personal thoughts).
Problem The issue is that you are using a continuous distribution 1,2,3 etc... for catagorical data and will result in confusing homozygotes with heterozygotes.
For example

If you have allele 3 (father) and allele 1 (mother) the "progeny" wil
be "allele" 2
If you have allele 2 (father) and allele 2 (mother) the "progeny"
will be "allele" 2

However, one is homozygous and the other is heterozygous and the allelic background for each is quite differnt.
The alternative is to use letters rather than numbers. You can always use "OneHotKeyEncoder" to convert back from AA, aa etc... to numerical catagories for analysis later (its a standard feature in Python).
If you used letters rather and then converted to numbers every "progeny" would have its own numerical catagory (but you would lose homozygous vs. heterozygous classification). It is standard practice to switch from letters to numbers and the "encoder" will ensure no mistakes are made. The other way around this is to designate the variable as a string, rather than as an integer... this might be more to your liking, but you are stuck in base 10.
With regards the throughput, 100000s 1000 individuals, a 100 million SNPs is doable via R, the operations are simple.
Just to clarify: it would be easy to script to a solution (even from base code), but I'm not sure about how to interpret the output. All I'm suggesting is you have another think about the goal of this analysis.

Again, my lack of clarity - 0,1,2 relates to B allele counts of bivariate SNPs (A/B) - these are not categorical values

In that case you appear to be violating Mendelian inheritance. Children can't receive 0.5 allele, which is what you are implying

No - they either receive 1, or 0. But on average, they receive 0.5.


text:
Question about Co-expression analysis and finding targets for lncRNAs
I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)
I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.
Few doubts:
1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?
2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?
3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.
Thank you

I think you should ask a separate question for your third question (remember to show what you have tried/searched...). For the question 1: which is your biological question? Are you interested in correlation between lncRNAs and mRNAs regardless if there is a tumor or not? Or are you interested in knowing if the co-expression in normal samples are different from the tumor? For question 2: You should be aware that you can find lots of false positive on those correlations. You should at least restrict yourself to those significantly correlated. But I wouldn't recommend to follow this path...

@llrs I'm interested in co-expression in tumor, may be I will do spearately. For the second, ofcourse I go based on pvalue.

In this case it is better if you only use the tumor samples. If you go based on p-value consider correcting the p-value to reduce the amount of false positives

And for WGCNA, after module detection, there are some important genes Im interested in were kept in the grey module which is unassigned genes. What to do in this situation?

You can then say that these genes are not co-expressed following a clear pattern. Either they don't follow a pattern, or are the ones regulating the expression of the other modules, or you don't have enough power/samples or it is too noisy to find a signal... The important genes might be for other reasons not for co-expression between mRNA and lncRNA (or they still are but you can't detect it)

You have several options to approach this with WGCNA (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.
Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).
I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.
I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.

Thanks a lot for the answer Peter. With WGCNA, I detected the modules and from that modules having mRNAs and lncRNAs how to detect hub genes inside modules?

Hub genes are the ones with high intramodular or eigengene based connectivity. Use the function signedKME to calculate kME values (it's really just a correlation of module eigengenes with individual expression profiles) and choose the top genes.

There is no hard and fast criterion for what constitutes a hub gene; I'd probably look at the first maybe 20 genes for a smaller module (less than say 200 genes) and maybe a few more (say up to 50) for larger modules; if you have a module of more than a thousand gene, looking at the top 100 may be useful. You can also go by kME (i.e., correlation); anything higher than 0.7 is pretty strong, and sometimes you can go down to say 0.6 to still find useful answers.

Thanks peter. After module detection, I used this code. Do you think this is right?         MEList = moduleEigengenes(datExpr, colors = dynamicColors);
MEs = MEList$eigengenes;
datKME=signedKME(datExpr, MEs);
head(datKME);
FilterGenes= abs(datKME$kMEbrown)>.8;

table(FilterGenes);
dimnames(data.frame(datExpr))[[2]][FilterGenes];

Hi Peter, Could you please check the above code whether that is right way to get hub genes?

As far as I can tell (the formatting is messed up) it should correctly give you the top hubs in the brown module.

A cursory search to bedtools documentation will reveal the bedtools closest feature - which might be exactly what you are looking in your third question. You can download the list of protein-coding genes from Ensembl biomart or UCSC table browser in bed format, and also convert the list of your lncRNAs into a bed compatible file. Rest is easy:
bedtools closest -a lncRNAs.bed -b proteinCoding.bed [-iu/-id] -D ref

2) If not WGCNA, if I use cor function in R with method Pearson, on
  what cutoff I should select the target genes?

Good call, this is a nice approach. It is a good question because the central issue is that a purists approach of 0.7 results in insufficient information to understand the behaviour of the data.
The guys here were using 0.4, in my experience co-variance of 0.5 and above is good evidence, whilst 0.7+ is strong evidence.
I've done this (not on RNAseq) and I used 0.25 to 0.45 and 0.45 and above then used a scatterplot with one half of the graph showing low co-variance and the other half strong co-variance. So in one sense I ignored my own advice (0.5+). In my defence its easy to reassign the boundaries. Each scatterplot (1000s) was then subject to tSNE analysis to assess the group behaviour. Its a different context but gives an idea.
Networks are cool for visualisation but they are not particularly  quantitive for most (not all) network theory.

Just to explain to everyone else, what happens if you get boundary effects with the data and you are assigning a discrete boundary on a continous distribution (its a negative binomial distribution in my case). I observed a vast amount of results  occurring between 0.2 to 0.25 - but 0.2 isn't good evidence, so that data has to be lost. These sort of boundaries were in days before genomic throughput. Thus in my context I didn't consider a hard boundary particularly informative.


text:
RSeQC provided bedfile does not work with read_distribution.py
I frequently make use of RSeQC, I used to use it with my own bed file, generated using gtf2bed, but now I thought I'd use the one provided by RSeQC themselves from here, I got the file hg38_Gencode_V28.bed.gz, when looking at it with IGV it looks similar to my self generated bedfile however, read_distribution.py (a script in RSeQC), does not report any reads in any of the categories (CDS_exons, Introns, etc).
Any idea why?
Edit: I see that my own bed file (from the GTF) had ENST identifiers like this ENST00000466557 whereas the file from RSeQC supplies them like ENST00000611746.1 (with decimals) that seems to be the only difference as far as I can see.

Ok, I figured it out already, it's and old issue I come accross every now and then. The original bed file lists the chromosomes as chr1,chr2,...,chrX whereas my bam files uses 1,2,...,X
I fixed the bed file like this (python3):
import pandas as pd
bedFile = '/hg38_Gencode_V28.bed'
bedFile = pd.read_csv(bedFile, sep='\t', header=None)
bedFile[0] = [value.strip('chr') for value in bedFile[0]]
newBedFile = 'hg38_Gencode_V28_reformatted.bed'
bedFile.to_csv(newBedFile, sep='\t', header=None, index=None)

And now read_distribution.py works.

You can also use chromosome mappings provided by @dpryan72 here: https://github.com/dpryan79/ChromosomeMappings


text:
Selecting 65000 SNPs where AF is close to 0.5 in all or most populations
I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.
My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.
I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.
My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:
gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /\;AF=0\.[34567]/) { if ($_ =~ /\w+\t\d+\t\S+\t\w\t\w\t/) {print $_} }'

Any recommendations?

If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.

I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.

What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."

My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).

Just use bcftools view for filtering:
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf

To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.
$ bcftools view -h input.vcf.gz|wc -l
255

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf

Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:
zgrep -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

And to limit to the first 65000 results:
zgrep -m 65000 -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

Finally, to keep the header of the original vcf unchanged, you could do:
zgrep -m 65000 -P '^(#|[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7])' my.vcf.gz

I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.
If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.

I have a different opinion about that. hg19 is outdated sind 2013. Whenever you start something new, you should move to the current reference genome, which is by now hg38. There are lifted over versions of 1000 Genomes, ExAC and gnomAD which provide at least the frequency data by ensembl. As long as the frequency data are given in the INFO column you can easily filter by them using bcftools. So there is no need to script something.

@finswimmer not sure what you are referring to, but if it's about hg38, I am not saying hg19 is better (it isn't) only that the gnomAD data is based on hg19 and liftover is imperfect. As for the rest, yes of course bcftools is the way to go!