text:
how to merge two fasta based on full fasta header and remove the duplicate
I am wondering if anyone can help to merge two fasta file based on the full header as initially gene ID may be similar in some of the sequences and to remove the duplicate entry based on full fasta header
I do have two multifasta file:
file1:
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

file2:
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

Expected output:
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

Please [edit] your question and show us i) an example of your input files including "matching" and "non-matching" sequences and ii) the output you would expect from those input files.

Thanks for the edit! Will sequences with the same header be identical? If so, why not extract all unique sequences directly instead of trying to find identical headers?

yes, same header have identical sequences, it will be great to extract unique sequences.

You can use the FastaToTbl and TblToFasta scripts I have posted here for this:
FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa

If you run this on your example, you will get:
$ FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa
$ cat merged.fa
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATG
GTGCGGGAACGCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGT
TTGGATGTCTCTGGCGACGG
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGC
GAACGCTACGGGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGAT
GTCACCGACGACGGTGGGTG
>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAG
GTCGCGGTCGGCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTC
GACCGCCTCGAAGCCCGCCT
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAG
GTCGCGGTCGGCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTC
GACCGCCTCGAAGCCAGCGC
>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGC
CAGTTGGTACGGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGAC
GTCGGCCTCGATGTCACGGG

This will remove any cases where both the header and the sequence are identical. If you instead want to do this based on the header only, not taking the sequences into account, you can just process both files and only print sequences from the 2nd if their ID line hasn't been seen before. For example, with awk:
$ awk '{
        if(NR==FNR){
            if(/^>/){seen[$0]++;}
            print
        }
        else{
            if(/^>/){ skip=seen[$0] }
            if(!skip){ print }
        }
       }' file1.fa file2.fa
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

it will be great to remove identical sequences as well based nucleotide sequence. but it is showing error: $ chmod +x FastaToTbl
$ chmod +x TblToFasta
$FastaToTbl file1.fasta file2.fasta | sort | uniq | TblToFasta > merged.fa
bash: FastaToTbl: command not found...
bash: TblToFasta: command not found...

@Dr.Yogesh yes, that's because you didn't put the scripts in a directory that's in your PATH. If you have them in the current directory, use ./FastaToTbl file1.fa file2.fa | sort | uniq | ./TblToFasta > merged.fa

Thanks, I am trying to copy the script again!

Thanks, it is working perfectly now, it was sensitive to extra space!


text:
Strange p-value histogram for differential gene expression analysis
I'm trying to perform pretty standard differential expression analysis using RNA-seq data.
I've used Kallisto to perform RNA quantification and am using Sleuth to perform the differential expression analysis.
I have 2 conditions that I'd like to compare, which is stored in the metadata:
> metadata
    sample cell replicate batch                   path
1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
6:  Test-3 Test         3   3   /path/to/Batch3/Test-3

I'm attempting to test for differentially expressed genes between the Base and Test cells while controlling for batch (the sequencing batch). The R code to do this is adapted from Sleuth's manual:
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
    aggregation_column = "ens_gene"
)

# smooth raw Kallisto abundances
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# fitting error measurements on reduced model
sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")

# compare models to calculate differentially expressed transcripts
sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")

# extract results
sleuth_table = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE
))
sleuth_table_tx = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE,
    pval_aggregate = FALSE
))

To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:

Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?
Edit following from @swbarnes2's answer:
I've used fdrtool to create an empirical null model of Sleuth's likelihood ratio test.
library(fdrtool)
local_fdr_txs = fdrtool(so$tests$lrt$`reduced:full`$test_stat, statistic="normal")

Which gives the following plots:

Plotting the distribution of these p-values produces:

This still has a hump, but looks more anti-conservative than before.
My remaining question is about whether I can use fdrtool to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test

what do mean by null hypothesis as normal ? Anyway, a good question +1.

I found that using the Wald test instead of the Likelihood Ratio Test produced very different results with respect to the p-values.
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    read_bootstrap_tpm = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes,
    aggregation_column = "ens_gene"
)

# same model of covariates
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# using Wald test instead of LRT
sleuth_obj = sleuth_wt(sleuth_obj, "cellLT1", "full")

# extract results
sleuth_table_genes <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = TRUE))
sleuth_table_txs <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = FALSE))

The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.

Try the protocol here for generating a different (smaller) dispersion estimate
https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues

I found this previously, but where DESeq2 has its test statistic be normal, Sleuth detects differentially expressed transcripts using a likelihood ratio test. fdrtool can create an empirical null model on normal data, but is the test statistic from Sleuth's test normal?

source: Stack Exchange

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

sort: votes


text:
Is there a standard tool used to convert a VCF to a BEDPE?
Many popular SV callers output a VCF. Unfortunately, there isn't a unified system at the present to label events with the same notation.
However, is there a standard method for converting these VCFs to BEDPEs?
svtools comes to mind: https://github.com/hall-lab/svtools
There are also some Bioconductor packages.
However, is there currently a standard widely used? Normally this is because a certain consortium did it in a certain way.
EDIT: There was a question in the comments about VCF vs BEDPE for structural variants.
https://github.com/jaquol/blog-like-old/blob/master/01-How-interpret-structural-variants-BEDPE-format.md
This is a structural variant in VCF format:
##fileformat=VCFv4.2
##source=LUMPY
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=STRANDS,Number=.,Type=String,Description="Strand orientation of the adjacency in BEDPE format (DEL:+-, DUP:-+, INV:++/--)">
##INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description="Imprecise structural variation">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">
##INFO=<ID=CIPOS95,Number=2,Type=Integer,Description="Confidence interval (95%) around POS for imprecise variants">
##INFO=<ID=CIEND95,Number=2,Type=Integer,Description="Confidence interval (95%) around END for imprecise variants">
##INFO=<ID=MATEID,Number=.,Type=String,Description="ID of mate breakends">
##INFO=<ID=EVENT,Number=1,Type=String,Description="ID of event associated to breakend">
##INFO=<ID=SECONDARY,Number=0,Type=Flag,Description="Secondary breakend in a multi-line variants">
##INFO=<ID=SU,Number=.,Type=Integer,Description="Number of pieces of evidence supporting the variant across all samples">
##INFO=<ID=PE,Number=.,Type=Integer,Description="Number of paired-end reads supporting the variant across all samples">
##INFO=<ID=SR,Number=.,Type=Integer,Description="Number of split reads supporting the variant across all samples">
##INFO=<ID=EV,Number=.,Type=String,Description="Type of LUMPY evidence contributing to the variant call">
##INFO=<ID=PRPOS,Number=.,Type=String,Description="LUMPY probability curve of the POS breakend">
##INFO=<ID=PREND,Number=.,Type=String,Description="LUMPY probability curve of the END breakend">
##ALT=<ID=DEL,Description="Deletion">
##ALT=<ID=DUP,Description="Duplication">
##ALT=<ID=INV,Description="Inversion">
##ALT=<ID=DUP:TANDEM,Description="Tandem duplication">
##ALT=<ID=INS,Description="Insertion of novel sequence">
##ALT=<ID=CNV,Description="Copy number variable region">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=SU,Number=1,Type=Integer,Description="Number of pieces of evidence supporting the variant">
##FORMAT=<ID=PE,Number=1,Type=Integer,Description="Number of paired-end reads supporting the variant">
##FORMAT=<ID=SR,Number=1,Type=Integer,Description="Number of split reads supporting the variant">
##FORMAT=<ID=BD,Number=1,Type=Integer,Description="Amount of BED evidence supporting the variant">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  3128_T47D
chr1    1584528 11  N   <DUP>   .   .   SVTYPE=DUP;STRANDS=-+:6;SVLEN=66541;END=1651069;CIPOS=-449,10;CIEND=-10,503;CIPOS95=-374,0;CIEND95=0,366;IMPRECISE;SU=6;PE=6;SR=0;PRPOS=2.44491e-07,3.26458e-07,4.33697e-07,5.7608e-07,7.57136e-07,9.87322e-07,1.28796e-06,1.66553e-06,2.14862e-06,2.75607e-06,3.52645e-06,4.48654e-06,5.68684e-06,7.14253e-06,8.95608e-06,1.11311e-05,1.38146e-05,1.70081e-05,2.09074e-05,2.54614e-05,3.0947e-05,3.7325e-05,4.47552e-05,5.33515e-05,6.32491e-05,7.45332e-05,8.7438e-05,0.00010173,0.000117821,0.000135629,0.000155258,0.000177125,0.000200949,0.000227194,0.000255827,0.000286183,0.000319177,0.000353501,0.000390203,0.000428746,0.000468895,0.000511378,0.000555476,0.000600756,0.000648001,0.000695831,0.000745204,0.000795631,0.000845517,0.000896997,0.00094752,0.000998877,0.00104988,0.00110053,0.00115056,0.00120038,0.00124937,0.00129849,0.00134617,0.00139206,0.00143762,0.0014814,0.00152386,0.00156521,0.00160519,0.00164404,0.00168131,0.00171717,0.00175139,0.00178475,0.00181657,0.00184663,0.00187563,0.00190327,0.00192975,0.00195479,0.00197842,0.00200183,0.00202369,0.00204549,0.00206583,0.00208582,0.00210445,0.00212224,0.00213857,0.00215452,0.00216896,0.00218263,0.00219584,0.00220845,0.00221975,0.00223061,0.00224049,0.00224956,0.00225841,0.00226617,0.00227317,0.00227931,0.00228505,0.00228996,0.00229495,0.00229948,0.00230347,0.00230755,0.00231109,0.0023146,0.00231767,0.00232087,0.0023235,0.00232593,0.00232843,0.00233037,0.00233266,0.00233471,0.0023365,0.00233822,0.00233974,0.00234121,0.00234261,0.00234391,0.00234524,0.00234656,0.00234787,0.00234901,0.00235006,0.0023513,0.0023522,0.00235324,0.00235429,0.00235519,0.00235614,0.00235695,0.00235786,0.00235864,0.00235957,0.00236033,0.00236117,0.00236198,0.00236288,0.00236365,0.00236441,0.00236514,0.00236581,0.00236671,0.0023675,0.00236851,0.0023693,0.00237023,0.00237106,0.00237196,0.00237283,0.00237371,0.00237471,0.00237556,0.00237647,0.00237733,0.00237819,0.00237897,0.00237973,0.00238056,0.00238145,0.00238222,0.00238294,0.00238394,0.0023848,0.00238545,0.00238618,0.0023871,0.00238796,0.0023889,0.00238976,0.00239053,0.00239142,0.00239222,0.0023932,0.00239412,0.00239504,0.00239599,0.00239704,0.00239815,0.00239902,0.00239984,0.00240073,0.00240184,0.00240288,0.00240385,0.00240489,0.00240597,0.00240695,0.00240787,0.00240872,0.00240953,0.00241046,0.00241129,0.00241228,0.00241313,0.00241411,0.00241502,0.00241603,0.00241712,0.00241811,0.00241934,0.00242034,0.00242141,0.00242256,0.00242369,0.00242488,0.00242598,0.00242723,0.00242845,0.00242958,0.00243082,0.00243182,0.00243282,0.00243389,0.00243485,0.00243593,0.00243699,0.00243818,0.00243921,0.00244038,0.00244162,0.00244278,0.00244386,0.002445,0.00244621,0.00244734,0.00244858,0.00244991,0.0024512,0.00245246,0.00245379,0.00245524,0.00245665,0.00245774,0.00245913,0.00246048,0.00246179,0.00246295,0.00246419,0.00246564,0.00246697,0.00246825,0.00246954,0.0024709,0.00247228,0.00247353,0.002475,0.00247639,0.00247812,0.0024795,0.00248075,0.00248215,0.00248359,0.00248517,0.00248655,0.00248793,0.00248927,0.00249072,0.00249223,0.00249351,0.00249512,0.00249648,0.00249801,0.00249976,0.00250133,0.00250293,0.00250464,0.00250618,0.00250758,0.00250932,0.00251113,0.00251256,0.00251406,0.00251567,0.00251749,0.00251919,0.00252051,0.00252197,0.00252351,0.00252518,0.00252675,0.00252834,0.00252988,0.00253136,0.00253307,0.00253469,0.00253652,0.0025383,0.0025403,0.00254212,0.00254387,0.00254575,0.00254754,0.00254942,0.0025514,0.00255316,0.00255534,0.0025573,0.00255935,0.00256123,0.00256313,0.00256486,0.00256668,0.00256838,0.00257033,0.00257235,0.00257408,0.00257612,0.00257812,0.00258004,0.00258208,0.00258394,0.00258613,0.00258805,0.00259035,0.00259238,0.00259449,0.00259664,0.00259868,0.00260091,0.0026033,0.00260538,0.00260742,0.00260951,0.00261187,0.00261382,0.00261579,0.00261808,0.00262029,0.00262253,0.00262481,0.00262695,0.00262929,0.00263173,0.00263396,0.0026362,0.00263851,0.00264069,0.00264344,0.0026458,0.00264821,0.00265024,0.00265247,0.00265488,0.00265735,0.00265965,0.00266185,0.00266402,0.00266663,0.00266919,0.00267144,0.00267383,0.00267621,0.002679,0.00268127,0.00268393,0.00268632,0.00268891,0.00269145,0.00269389,0.00269661,0.00269912,0.00270165,0.00270404,0.0027067,0.00270939,0.00271195,0.00271455,0.00271722,0.00271969,0.00272252,0.0027254,0.002728,0.00273082,0.00273369,0.00273648,0.00273952,0.00274227,0.00274511,0.00274816,0.00275103,0.00275394,0.00275687,0.00275977,0.00276241,0.00276495,0.00276767,0.00277064,0.00277351,0.00277641,0.00277909,0.00278202,0.00278452,0.00278729,0.00278998,0.00279275,0.0027958,0.00279878,0.0028016,0.00280472,0.00280771,0.0028109,0.00281415,0.00281717,0.00282038,0.00282346,0.00282626,0.00282919,0.00283207,0.00283518,0.00283843,0.00284122,0.0028439,0.00284708,0.00284994,0.00285294,0.00285587,0.00285848,0.00286126,0.0028643,0.00286701,0.00287003,0.00287284,0.00287589,0.00287886,0.00288176,0.00288469,0.00288762,0.0028909,0.00289392,0.00289664,0.00289945,0.00290244,0.00290543,0.00290824,0.00291089,0.00291354,0.00291634,0.00291902,0.00292172,0.00292413,0.0029271,0.00292984,0.002933,0.0029355,0.00293802,0.0029405,0.00294333,0.00294613,0.00117372,0.000467584,0.000186264,7.42093e-05,2.95642e-05,1.17781e-05,4.69226e-06,1.86919e-06,7.44696e-07,2.96676e-07;PREND=3.02026e-15,4.78471e-14,7.57973e-13,1.20067e-11,1.90178e-10,3.01282e-09,4.77324e-08,7.56151e-07,1.19779e-05,0.000189743,0.00300582,0.0030039,0.00300127,0.00299837,0.00299541,0.00299255,0.00299047,0.00298793,0.00298555,0.00298267,0.00298033,0.0029782,0.00297579,0.00297295,0.00297017,0.00296769,0.00296562,0.0029632,0.00296072,0.00295816,0.00295557,0.00295308,0.00295062,0.0029479,0.0029447,0.00294211,0.00293963,0.00293734,0.00293449,0.00293173,0.00292904,0.0029266,0.00292425,0.00292122,0.00291848,0.00291567,0.00291333,0.00291077,0.00290832,0.00290537,0.00290246,0.00289957,0.00289704,0.00289385,0.00289114,0.0028884,0.00288549,0.00288284,0.00287995,0.00287726,0.00287463,0.00287144,0.00286885,0.00286569,0.0028629,0.00286073,0.00285801,0.00285472,0.00285196,0.00284879,0.00284594,0.00284346,0.00284056,0.00283777,0.00283502,0.0028323,0.00282939,0.00282645,0.00282353,0.00282053,0.00281778,0.00281457,0.00281189,0.00280891,0.00280615,0.00280355,0.00280057,0.00279736,0.00279468,0.00279185,0.00278916,0.00278699,0.0027844,0.00278164,0.00277904,0.00277575,0.00277311,0.00277044,0.00276776,0.00276461,0.00276186,0.00275906,0.00275677,0.00275435,0.00275168,0.00274942,0.00274617,0.00274363,0.00274107,0.00273819,0.0027355,0.00273297,0.00273059,0.00272816,0.00272563,0.00272344,0.00272085,0.00271833,0.00271609,0.00271369,0.00271133,0.00270894,0.0027064,0.00270412,0.00270171,0.0026994,0.00269703,0.00269418,0.00269178,0.00268936,0.00268728,0.00268506,0.00268312,0.00268036,0.00267787,0.00267564,0.00267347,0.00267118,0.00266877,0.00266653,0.00266441,0.00266271,0.00266054,0.00265825,0.0026559,0.00265383,0.00265161,0.00264968,0.00264732,0.00264539,0.00264295,0.00264098,0.00263869,0.00263636,0.00263434,0.00263215,0.00263023,0.00262822,0.00262628,0.00262422,0.00262236,0.00262036,0.00261831,0.00261628,0.00261451,0.00261257,0.00261059,0.00260871,0.00260687,0.00260501,0.00260267,0.00260094,0.00259868,0.00259649,0.00259478,0.00259308,0.00259143,0.00258955,0.00258721,0.00258551,0.00258347,0.00258174,0.00258005,0.00257832,0.00257656,0.00257497,0.00257335,0.00257165,0.00256984,0.00256803,0.00256624,0.00256436,0.00256254,0.0025607,0.00255925,0.00255792,0.00255616,0.00255439,0.0025528,0.00255149,0.00254971,0.00254832,0.00254703,0.00254555,0.00254428,0.00254283,0.00254166,0.00254023,0.00253862,0.00253728,0.002536,0.00253436,0.00253252,0.00253096,0.00252943,0.00252819,0.00252656,0.00252508,0.00252379,0.00252243,0.00252084,0.00251975,0.00251831,0.00251666,0.00251521,0.00251391,0.00251267,0.00251144,0.00251024,0.00250897,0.00250758,0.00250646,0.00250525,0.00250396,0.00250248,0.00250084,0.00249962,0.00249836,0.00249712,0.00249567,0.00249435,0.00249341,0.00249217,0.00249093,0.00248988,0.00248849,0.00248716,0.00248589,0.00248488,0.00248398,0.00248264,0.00248155,0.00248047,0.00247944,0.00247824,0.00247724,0.00247609,0.00247492,0.00247375,0.00247262,0.00247147,0.0024706,0.00246888,0.0024678,0.00246675,0.00246572,0.00246487,0.00246414,0.00246316,0.0024619,0.00246083,0.00245995,0.00245894,0.002458,0.0024569,0.00245599,0.00245492,0.0024538,0.00245275,0.0024515,0.00245035,0.0024494,0.00244835,0.00244731,0.00244625,0.00244542,0.00244473,0.0024439,0.00244295,0.00244172,0.00244076,0.00243982,0.00243907,0.00243812,0.00243725,0.00243639,0.00243555,0.00243464,0.0024338,0.00243277,0.00243145,0.00243048,0.00242965,0.0024285,0.00242736,0.00242654,0.00242575,0.00242483,0.00242391,0.00242298,0.00242193,0.00242078,0.00241981,0.00241862,0.00241771,0.00241663,0.0024155,0.0024146,0.00241367,0.0024125,0.00241116,0.00241008,0.00240902,0.00240729,0.00240548,0.00240419,0.00240253,0.00240081,0.00239908,0.00239741,0.00239517,0.0023929,0.00239052,0.00238798,0.00238505,0.00238189,0.00237904,0.00237565,0.00237201,0.00236829,0.00236374,0.00235923,0.00235403,0.00234858,0.00234288,0.00233677,0.00233038,0.00232438,0.0023171,0.00230981,0.00230198,0.00229329,0.00228476,0.00227551,0.00226594,0.00225598,0.00224564,0.00223392,0.00222138,0.00220881,0.00219607,0.00218362,0.00217038,0.00215639,0.00214182,0.00212751,0.00211102,0.00209334,0.00207535,0.00205704,0.00203711,0.00201595,0.00199504,0.00197257,0.00194934,0.0019241,0.00189856,0.0018717,0.00184359,0.00181463,0.00178446,0.00175283,0.0017206,0.0016859,0.00165096,0.00161563,0.00157919,0.00154094,0.00150137,0.00146059,0.00141945,0.00137825,0.0013367,0.00129363,0.00124908,0.00120422,0.00115974,0.00111459,0.00106967,0.00102485,0.000979928,0.000935386,0.000891343,0.000847616,0.000805325,0.000763099,0.000721605,0.000680876,0.000641258,0.000602415,0.000564516,0.000529212,0.000494813,0.000461726,0.000429703,0.000398541,0.000369522,0.000342268,0.000316742,0.000292816,0.000270237,0.000249087,0.000228998,0.000210083,0.000192775,0.000176194,0.000160682,0.000146541,0.000133475,0.000121171,0.000110027,9.98983e-05,9.04025e-05,8.1898e-05,7.37485e-05,6.64247e-05,5.991e-05,5.34662e-05,4.79254e-05,4.28124e-05,3.80863e-05,3.37871e-05,2.98488e-05,2.64158e-05,2.31778e-05,2.02341e-05,1.77006e-05,1.53314e-05,1.32501e-05,1.13287e-05,9.72139e-06,8.17381e-06,6.84904e-06,5.73097e-06,4.74769e-06,3.91083e-06,3.20076e-06,2.6024e-06,2.08961e-06,1.69223e-06,1.36759e-06,1.08247e-06,8.47588e-07,6.69622e-07,5.17384e-07,4.00954e-07,3.11072e-07,2.34849e-07,1.76165e-07,1.31673e-07,9.73364e-08,7.19841e-08,5.32892e-08,3.90585e-08,2.83402e-08,2.02418e-08,1.4579e-08,1.06446e-08,7.59531e-09,5.37474e-09,3.79538e-09,2.64724e-09,1.84573e-09,1.29566e-09,9.10203e-10,6.14793e-10,4.26199e-10,2.87287e-10,1.936e-10,1.30732e-10,9.00388e-11,6.2542e-11,4.32799e-11,2.90205e-11,1.95416e-11,1.3365e-11,8.69465e-12,5.75237e-12,3.71099e-12,2.54394e-12,1.73591e-12,1.14175e-12,7.46542e-13,4.68409e-13,3.08728e-13,1.94589e-13,1.17053e-13,7.58227e-14,4.59278e-14,2.89681e-14,1.97988e-14,1.19182e-14,7.29947e-15,4.37438e-15  GT:SU:PE:SR ./.:6:6:0

And here is the same event as a BEDPE:
chr1    1584078 1584538 chr1    1651058 1651572 11  .   -   +   DUP .   SVTYPE=DUP;SVLEN=66541;END=1651069;STRANDS=-+:6;IMPRECISE;CIPOS=-449,10;CIEND=-10...

https://github.com/cancerit/BRASS/wiki/BEDPE
The following is a table describing how SVs are represented in BEDPE format:
Heading Type    Description
chr1    String  Chromosome of lower coordinate
start1  0-based int Start coordinate of lower coordinate
end1    1-based int End coordinate of lower coordinate
chr2    String  Chromosome of high coordinate
start2  0-based int Start coordinate of high coordinate
end2    1-based int End coordinate of high coordinate
id/name String  ID of event, correlates with VCF
brass_score int Number of aberrant pairs contributing to the rearrangement group.
strand1 [+-]    Strand of end in 'genomic' context - see table
strand2 [+-]    Strand of end in 'genomic' context

Do you really mean bedpe? I may be missing something, but what would converting a vcf to bedpe even mean? Bedpe is normally used to described paired sequences, usually primers. What would converting a VCF to bedpe entail? Could you maybe [edit] your question and give an example of a minimal VCF and what it would look like as bedpe?

@terdon I've edited the question for clarity

Thanks! That makes it much clearer. I don't really see how you'd represent deletions, but if the LUMPY folks are doing it, I guess it must be possible. Did you try vftobedpe from svtools?

SURVIVOR has a subcommand to convert an SV vcf to bed, although I'm not sure if it's actual bedpe format.
To install:
git clone https://github.com/fritzsedlazeck/SURVIVOR.git
cd SURVIVOR/Debug
make

Then:
./SURVIVOR vcftobed input.vcf $minSize $maxSize output.bed

If you run this on your example file, you get:
$ ./SURVIVOR vcftobed foo.vcf 1 10000000000000000 bar.bed; cat bar.bed
chr1    1584528 1584528 chr1    1651069 1651069 11  ,   -   +   DUP

Could you edit your answer and add this command's name and how to execute it? As it stands, this is kind of more of a comment than an answer.

Feel free to move it to a comment then. I consider the documentation of the tool sufficient and clear, and do not believe that copy-pasting that usage information is beneficial, as it might get outdated while the official documentation shouldn't.

It's a question of what is considered an answer on the SE platform. This seems like a good answer, but we really want the relevant information to be here in case the link goes stale and to be able to have the relevant details here. You're linking to the docs, so if they change and your instructions become obsolete, the docs will still be there. In the meantime, however, this is not actually giving us enough info to know how to use it and that's a shame.

Hi @terdon, thanks for the edits. "However, this isn't a bedpe file, so it looks like it isn't what the OP was asking for". That looks pretty close to a BEDPE, right? What is the difference you see, besides the file extension?

Whoops, my bad! It is indeed a bedpe and not a bed. I had managed to miss the second set of coordinates, sorry! However, it doesn't seem to be right for what the OP wants since that's just repeating the start position two times and then the end position another two times. Have a look at the blog post the OP linked to for what they seem to be looking for.

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


text:
Reducing a string to a certain amount of the first characters in R
So I have a list of several thousand RNA transcripts, and I am trying to find out how they are spliced. Right now I am running matchLRpatterns() from the Biostrings package with a max gap length of 0, after running a matchpattern function to categorize the transcripts by donor sites (where the first cut in an RNA transcript is made to cut out introns). But applying that to several thousand transcripts is quite time consuming, when you have 5 donor sequences to match to 9 acceptor sequences in 100 nucleotide long transcripts. Instead, I am thinking to sort by donor, remove all that comes before the splice junction by converting it to a string then using gsub, then using matchpattern again, but I only want to check the first 8 nucleotides after the splice junction for the consensus sequence that follows my list of acceptor sites. If I just run matchpattern without shortening the length of the transcript, I will just find acceptor sites in the RNA transcript, not the acceptor site at the splice junction.
I was thinking lapply the gsub function to the transcripts to remove up to the splice junction, then lapply the strsplit function to make a list of lists of groups of 1 character strings, lapply the head function to get the first 8 single character strings in each list, lapply the paste function to make them one string again, lapply the DNAString function to make the sequences readable by matchpattern, then lapply the matchpattern function to check the acceptor sequence at the splice junction. Is there a better way? This feels a bit clunky and roundabout way of doing it. I know gsub can remove characters after a certain character or sequence of characters, but I am not sure what the first sequence in each transcript is at the splice junction. Can I still use gsub, but tell it to remove characters after a certain number of characters? Or better yet, is there a way to get matchpattern to only check the first x nucleotides in a sequence?

Please [edit] your question and include i) an example of your input data, ii) the output you would want from that data and iii) the code you are currently using so we don't reinvent the wheel.

I suggest you follow @terdon suggestion to increase the likelihood of receiving help. The task  you are referring to in the title can be easily achieved via the following: smallx<-substr(x, start, stop)

Is there any particular reason why you want to do this via pattern matching in R? This sounds like something that SuperTranscripts and Lace might be useful for:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1284-1

Finally, Lace will annotate each superTranscript with blocks either by using the graph structure itself or alternatively using a script called Mobius (packaged with Lace) which infers the annotation from spliced junctions discovered when mapping reads to the superTranscript.

I was unaware of the existence of Lace and SuperTranscripts. My mentor in the lab I just started in wanted me to write code in R to help organize his Illumina reads for his research. I can ask him if he has heard of it, and if he would consider using it.


text:
Raw Data frame manipulation in python
Using python 3 I need to process qPCR sequencing raw data outputs by searching for the first occurrence of a user defined string and then making a new data frame using all lines after that string. I am trying to find solutions in the pandas doc but so far unsuccessful.
This is a raw output .csv file that I need to process. (couldn't paste complete csv as exceeds character limit, this is lines 40-50 and am hoping this text is useful?). I need to tell pandas to create a new data frame that 1. starts at the line containg the first occurance of str("Sample Name") with that line as header and containing all lines following. And then 2., only including columns ("Sample Name"), ("Target Name"), ("CT").
Could someone please help me so that I can use python to help me analyze biological data?
Many thanks,
Luke
40,Quantification Cycle Method,Ct,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

41,Signal Smoothing On,true,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

42,Stage where Melt Analysis is performed,Stage3,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

43,Stage/ Cycle where Ct Analysis is performed,"Stage2, Step2",,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

44,User Name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

45,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

46,Well,Well Position,Omit,Sample Name,Target Name,Task,Reporter,Quencher,Quantity,Quantity Mean,SE,RQ,RQ Min,RQ Max,CT,Ct Mean,Ct SD,Delta Ct,Delta Ct Mean,Delta Ct SD,Delta Ct SE,Delta Delta Ct,Automatic Ct Threshold,Ct Threshold,Automatic Baseline,Baseline Start,Baseline End,Amp Status,Comments,Cq Conf,CQCONF,HIGHSD,OUTLIERRG,Tm1,Tm2,Tm3,Tm4

47,1,A1,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.357698440551758,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9588544573203085,N,Y,N,81.40960693359375,,,

48,2,A2,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,24.05980110168457,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,15,Amp,,0.9592687354496955,N,Y,N,81.40960693359375,,,

49,3,A3,False,WT1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.012556076049805,23.4766845703125,0.5336655378341675,,,,,,True,20959.612776965325,True,3,16,Amp,,0.9592714462250367,N,Y,N,81.40960693359375,,,

50,4,A4,False,fla11fla12-1,AtTubulin,UNKNOWN,SYBR,None,,,,,,,23.803699493408203,24.419523239135742,0.5669151544570923,,,,,,True,20959.612776965325,True,3,17,Amp,,0.9671570584141241,N,Y,N,81.40960693359375,,,

This is the code that I have so far:
import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

Please [edit] your question and add i) an example of your input data, ii) an explanation of the data frame you want to build from that data and iii) the code you have so far so we don't reinvent the wheel.

Don't use grep using subprocess to filter dataframes. That's crazy :-D
Follow a pandas tutorial: the filtering that you want to do on your dataframe can be done with pandas commands, no need to use grep. If you are using python, avoid mixing your code with shell commands.

Please don't give us photographs of your data! We need your actual data so we can test any solutions we come up with. Just save the excel file as csv or tsv or something and post the first few lines here. And paste the errors as well. You can use the formatting tools to make all this look like code.

What is that header about in your csv file? I suppose that a csv file should be just a comma-separated table. What is the first column? line numbers? Why do they start at 32? I understand that starting with python might be difficult, have you tried some tutorials?

Thanks for all resposes. @terdon that formatting link helped,

@Kamil S Jaron I have tried to explain my csv example and what I need to acheive a bit better, hope that's helpful?

data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

This doesn't make much sense ... you read in Excel to a dataframe, you then read out to csv - the dataframe doesn't have an index, but you include an index and then read it in again. This could do more harm than good.
Anyway your supervisor is right it is very easy to parse this data using pandas, although there will be more expertise doing this via R.
I don't want to peer through the gloom of a CSV output, but I would like to highlight how easy it is to parse this stuff into a dataframe. The dataframe 'data' is fine.

Parse You can simply remove the top rows of bumf data = data.iloc[6:] (I might have miss counted, plus see below),
Headers You will probably have to remove the "headers" row to a list and
then introduce back onto the dataframe using data.columns =
['header1', 'header2' ...] ... rather this is the list (array in Perl) which you can automatically generate using a magic to_list type pandas command which I forget, its not list - its cuter
Columns You next might use drop to remove the columns you don't want (that might be alot of columns) or else make a subset of
the dataframe into a new dataframe, you can also use iloc for
columns too.
Index You then assign A1, A2, ... etc  as an index data.index = ['A1',
'A2', ]...etc..  you can reshape column 3 to automate this
to_numeric You might have to assign the CT values as floats
    or better to_numeric
Extrapolation You would then define a subroutine (or whatever they call them in
Python) so perform the extrapolation and use the apply command on
your dataframe and assign it to a new column, where the input is your CT values.
Output Finally, issue your to_csv command to_pickle is also a good idea because you have a df that is worth importing at a future date.
Plotting this is easy using data.plot which appears to be some sort of inbuilt matplotlib, keep in mind seaborn DONT try using R's ggplot2 within Python using rpy2 .... its a NIGHTMARE to get working and doesn't support the latest ggplot2 (!!!). R
s ggplot2 is stronger  than seaborn or matplotlib (providing you are working within R), but these two are good for quick graphs.

If I had time I'd write the code for you, but I'm afraid I'll have to leave it there. I'm still in debt with a Perl efetch script I promised a previous OP.
In summary, Python pandas is fabulous stuff for this purpose, I love it. Doing this via Perl and "subprocessing" inside R... its doable, well in this particular case I suspect its just as easy, but generally its alot more code and often fiddly.

Really appreciate your help! Thanks. When I get a chance I will post my updated code.


text:
How can I correctly get 3D scatter plot for micro-array pheno data in R [updated]?
I intend to see data points distribution (each reows) within labeled groups (different batchs such as 1,2,3, and so on) in 3D scatter plot, because I want to see the distribution of the data points and want to see how similar each group of data points in 3D space. To do so, I used scatterplot3d package from CRAN to get 3D scatter plot, didn't get correct plot for my data.
reproducible data:
Here is the reproducible data that simulated from the original dataset.
persons_df <- data.frame(ID= LETTERS[1:20], ages=sample(1:50,20, replace = FALSE),
                    batch=c(rep(1,5), rep(2,6),rep(3,4),rep(4,5)),
                    platform=c(rep("CATH",12), rep("ANN",8)))

basically, I want to render 3D plot for the data points that belongs to the same batch with the same color (1 is red, 2 is blue, 3 is green, 4 is yellow).
this was my attempt to render PCA plot
require(scatterplot3d)
par(mar=c(4,4,4,4), cex=1.0, cex.main=0.8, cex.axis=0.8)
scatterplot3d(persons_df$batch)

but above code doesn't work because first I should group data points that belongs to different batches then render them in 3D space. This is not very intuitive to me how to do it in R? Any idea to render possible 3D scatter plot for my data? any thought?
desired 3D plot:
I want 3D plot something like this:

how can I make this happen? any idea to get this done?

PCA doesn't make sense here. PCA is a dimensional scaling technique. When we use PCA, we are trying to represent (usually) highly dimensional data in two or three dimesnions. For example, in the example you linked prcomp(iris[, -5]) is performing PCA on sepal length, sepal width, petal length, and petal width. So we have 4 dimensions here.
In your reproducible example, you only have two dimensions to look at, Age and Batch. Number of dimensions aside, I think you are missing the useful data on which it would be useful to apply PCA. This is just metadata.

@KohlKinning agreed, I am going to rephrase my question. I need 3D scatter plot to see data points (a.k.a, number of rows that belongs to different batch). Do you have any idea?


text:
Error: package or namespace load failed for 'ChIPseeker'
How do I solve the problem of loading the  bioconductor Chipseeker package? I have installed in the following way and I having an error while loading it and I do not know why.
BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

library(ChIPseeker) # produces and error

Error: package or namespace load failed for 'ChIPseeker':
  'TxDb.Hsapiens.UCSC.hg19.knownGene' is not exported from 'namespace: TxDb.Hsapiens.UCSC.hg19.knownGene'

Thank you in advance for any help on how to go about this.
session info:

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] BiocInstaller_1.32.1                    TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [3] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.0 GenomicFeatures_1.34.8
 [5] AnnotationDbi_1.44.0                    GenomicRanges_1.34.0
 [7] GenomeInfoDb_1.18.2                     IRanges_2.16.0
 [9] S4Vectors_0.20.1                        cellrangerRkit_1.1.0
[11] Rmisc_1.5                               plyr_1.8.4
[13] lattice_0.20-35                         bit64_0.9-7
[15] bit_1.1-14                              ggplot2_3.1.1
[17] RColorBrewer_1.1-2                      Biobase_2.42.0
[19] BiocGenerics_0.28.0                     Matrix_1.2-14

loaded via a namespace (and not attached):
Error in x[["Version"]] : subscript out of bounds
In addition: Warning message:
In FUN(X[[i]], ...) : DESCRIPTION von Paket 'digest' fehlt oder zerstört

There might be a problem with your current R session if sessionInfo() gives you an error. Try restarting the R session (if you're using RStudio) or just exiting and starting a new R session (if you're on the command line).

Simple, use the latest version of R and Bioconductor, now on version 3.6.1 and Bioconductor 3.9. Those are the versions that are supported. If you can't upgrade R, check that BiocManager::valid() is true, you might have ended in some inconsistency between Bioconductor versions.
Looking at your sessionInfo, you seem to have some problems with the installed libraries, check and correct those too if possible

@IIrs: BiocManager::valid() gives the Warning message: 10 packages out-of-date; 1 packages too new.  How can I know the package which is too new?

It also printed some more messages. Please edit the original question and post the full comment. If it didn't (I don't know how that could be) use BiocManager::install() (without any package), it should update the out-of-date packages. But the problem might be that this too-new is the one you have problems with. Store the result of a <- BiocManager::valid() and explore the a$toonew(or something alike).


text:
Error while calling bcftools mpileup - Failed to open -: unknown file type
I have sequenced a bacterial genome with a GridIon from ONT. Basically what I want to check is whether or not trimming 50 bps at the beginning of the reads will improve alignment against the reference genome and ultimately the call of a consensus sequence.
For that what I have done so far is:

concatenated all the fastqs generated from the Gridion into a single (very large) fastq file: ref494_cat_all.fastq
aligned this concatenated fastq against the reference genome ref_genbank.fasta with minimap2: ~/Downloads/minimap2-2.17_x64-linux/minimap2 -ax map-ont ref_genbank/ref_genbank.fasta ref494_cat_all.fastq > ref494_aligned_trim50.sam
Performed the general samtools steps to covert sam to bam, index and sort

Then, in order to get the consensus sequence I first try to generate a .vcf.gz file, according to this tutorial
bcftools mpileup -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz

Immediatelly after calling the command above I got:
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[mpileup] 1 samples in 1 input files

Then, after a minute or two, I got:
Failed to open -: unknown file type

I'm pretty sure that there's no issue with .bam file, as I could visualise the alignment with IGV or tablet.
Does anyone know what the issue can be?

Please capture the output of the first command and show us the first few lines (via bcftools view if it's not text): bcftools mpileup -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam > output.tmp

I got Killed output to the terminal after a minute or two. Also this command generates an empty output.tmp file.

As noted below, killed due to memory exhaustion. However: you deserve an error message. Mpileup typically doesn't use that much memory and I have been unable to reproduce this on (non-ONT) data here. If you are able to, it would be interesting to try current develop bcftools to see if this still happens and/or to run it under the debugger and see where it crashes.

The latter is quite easy: run gdb --args bcftools mpileup -f ref_genbank/ref_genbank.fasta -o /dev/null ref494_aligned.sorted.bam, then type run and wait for it to crash. Then type backtrace; the output from this will be instructive.

Just before the program got killed I got Program terminated with signal SIGKILL, Killed. The program no longer exists.
The output of backtrace was No stack
That you cannot reproduce this error on non-ONT data, you think it can be related to the size of my input files? My bam is 14GB (coverage is about 2.000, we ran for about 48h, which we won't do routinely)

Oh yeah, SIGKILL. Perhaps you were hit by the OOM Killer (cf https://stackoverflow.com/questions/12288550/c-linux-binary-terminated-with-signal-sigkill-why-loaded-in-gdb) so the lack of error message maybe isn't bcftools's fault. 14G isn't a big BAM file, but the level of coverage and nature of ONT data (long reads with maybe interesting CIGARs) would be things that could make mpileup take more than the usual amount of memory.

Just tried to run the original command (as in my question) on a cluster (24 cores, 48GB of RAM). It ran for an hour or so, then I got the same error Failed to open -: unknown file type, so I guess I will try some commands without piping to bcftools to see what output I will get

The - is usually used to mean standard input when reading data. So the message would suggest that the second bcftools command fails to read from stdin, so it fails to read the output of the first bcftools command:
$ echo "foo bar" | /ref/software/bcftools-1.8/bcftools call -mv -Oz -o calls.vcf.gz
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
Failed to open -: unknown file type

So there is something wrong with either your input file or with the output produced when running the first bcftools command on that input file. We can't really help more unless you share the file.

I doubt there's any problem with the input files, as I can visualise the alignment on table or IGV, so I guess the bcftoools mpileup is not generating any output (please see my comment to @John Marshall). I guess I will report the problem on their github page.

@BCArg yes, precisely. If you got Killed in the terminal, then the process is being killed by the operating system, most likely because it is taking up too much memory. Sounds like you need to run this on a more powerful machine.

yes mystery solved, I monitored the the performance of the command bcftools mpileup -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam > output.tmp with htop and the program got killed right after I reached the maximum of my RAM memory (16GB), so in fact running the command on the cluster appears to be the solution.

It's possible you're looking at instructions for a recent version of samtools / bcftools, whereas your computer has an older version. Can you try samtools mpileup -g -f ... to output the pileup as BCF format instead of bcftools mpileup -f ...?
Update: it's good to know that the samtools / bcftools pipelines are at least producing consistent results, which suggests that the error is not associated with the output format on the mpileup side of the pipeline, and also good to see that you've identified that this was a memory issue.
The most frequent cause of memory issues I've had with sorted BAM files has been trying to figure out variants in extremely high coverage data (e.g. over 1000X), which would make sense for a bacterial isolate sequenced on a GridION. I use my own script to carry out a reservoir sampling on mapped BAM format files, which reduces coverage down to a more manageable level. Here's an example command to normalise to 100X coverage where possible:
samtools view reads.bam | ./samNormalise.pl -coverage 100 | \
   samtools view -b > reads_normalised.bam

did not work unfortunately, got the following messages after calling:
samtools mpileup -g -f ref_genbank/ref_genbank.fasta ref494_aligned.sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz
[warning] samtools mpileup option g is functional, but deprecated. Please switch to using bcftools mpileup in future.
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[mpileup] 1 samples in 1 input files
Failed to open -: unknown file type

The VCF/BCF-outputting mpileup command was transferred from samtools to bcftools relatively recently. So in fact it is @gringer's suggestion that is instructions for an earlier version of samtools / bcftools.


text:
Is there a computational method/method/way to predict if a cell population in a tissue immigrated or is enriched locally?
I have single cell data which I have analysed for differential expression. In my experiment,  I subjected two groups of mice (a control and a treatment group) to treatment that will lead to the proliferation of my cells of interest in the treated mice. Other scientists have said that the increase of cell numbers caused by this treatment results from increased recruitment from precursor cells from the bone marrow while other scientists argue that the increase in cell number is due to local proliferation.
I wish to address this question computationally, based on my single cell rna seq data using R or any other method. I have googled but I can find no suitable tool. Does it suffice to test the expression of known cell migration markers, produce a migration score  and judge based on the migration scores of my proliferating cells?
If not, I will be very grateful for suggestion of a computational tool, method or package that can help me address this question even indirectly.
Thank you in advance.

Do you have single cell data of both regions, bone marrow and local region? How good are your cell migration markers?

@Irls Thanks a lot for your response. Yes I do have but do you please have a hint on now to proceed from here?

This may not be a viable answer for all experimental contexts, but you can compare expression profiles of individual cells against datasets that contain known reference cell types (migratory and not). We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations in scRNAseq experiments. There are various reference datasets available, but you can also provide a custom reference by preparing data from databases like GEO.
You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. If your reference contains both migratory and non-migratory cell types, maybe CIPR can help differentiate them. For other identity prediction tools, check out this article is a good resource that compares different cluster calling algorithms.


text:
What do the read colors in IGV mean?
I was looking at a bam file in the IGV viewer and saw:

What do the different read colors mean? Why is one read a light blue color, another green, another aquamarine (?), another purple and another two orange? What do these different read colors represent? The IGV docs have sections on interpreting color by insert size, pair orientation, bisulphide mode, splice junctions and sashimi plot, but none of them seems to be describing what I'm looking at there.
More generally, where are all these default behaviors documented? Are they somewhere in http://software.broadinstitute.org/software/igv/UserGuide and I'm missing them?

Could it be that these are reads with their mates found on other chromosomes? That's the part of the IGV documentation which seems to make sense in this context:

for paired end reads that are coded by the chromosome on which their mates can be found


text:
Why is 42 the seed used in RunPCA() in Seurat?
The RunPCA() in Seurat 3.0 function has seed.use=42 by default. How does one choose a seed number? Is this chosen simply because it is the Answer to the Ultimate Question of Life, the Universe, and Everything? Or maybe Lewis Carroll's curious obsession with the number?

Any integer >= 0 (but not too large) can make a good seed, 42 is a popular choice for the reasons you note, https://www.residentmar.io/2016/07/08/randomly-popular.html and https://blog.semicolonsoftware.de/the-most-popular-random-seeds/

The authors of the Mash and sourmash algorithms use 42 as the seed for the MurmurHash3 k-mer hashing function, and I understood this to be a nod to the Hitchhiker's Guide as you suggested. So the idea isn't completely without precedent.

Given that it is virtually impossible for a human to predict what random numbers will be generated given a certain seed, and what effects they will have for a given application, the choice of the seed can be considered arbitrary.
What matters is that fixing the seed should fix the behaviour of the application. What seed has been chosen does not really matter, as long as it is documented.
The choice of 42 could therefore well come from one of the reasons you propose. Only the authors of the programs could confirm.


text:
How are UMIs used to dedupulicate in Drop-seq tools?
The module DigitalExpression which is part of the popular Drop-seq tools digitally count gene transcripts.
The manual is not very clear on how exactly it resolves the UMI information.
Maybe I'm missing something but I was wondering if UMI are used in isolation or in conjunction with the barcode and or transcript?
The manual i.e. Drop-seq Alignment Cookbook only states:
To digitally count gene transcripts, a list of UMIs in each gene, within each cell, is assembled, and UMIs
within edit distance = 1 are merged together. The total number of unique UMI sequences is counted,
and this number is reported as the number of transcripts of that gene for a given cell.

To clarify my point, I'll use an example. Say my reads are only 4 bp and read 1 has the cell barcode at position 1+2 and UMI at position 3+4.  read 2 constitutes the transcript read.
Given the following example reads:
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = AATT  read2 = TTTT

UMI and barcode are the same but read2 is different, would they be collapsed to 1 UMI count or 2?
Similarly as case with same UMI and transcript but different cell barcode:
seqence1:   read1 = AATT  read2 = ATGC
seqence2:   read1 = GGTT  read2 = ATGC

Also of note, UMI deduplication is the default setting as of the parameter below:
OUTPUT_READS_INSTEAD=Boolean  Output number of reads instead of number of unique molecular barcodes.  Default value:
                              false. This option can be set to 'null' to clear the default value. Possible values:
                              {true, false}

UMI and barcode are the same but read2 is different, would they be
  collapsed to 1 UMI count or 2?

If the two reads share cell barcodes and UMI sequences and align to the same gene, (even if they are different sequences) then they originate from the same RNA molecule, and should only count as 1.
https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry

Do you have any sources (documentation, code, citation etc) to back this up? To my knowledge the UMI is designed to spot PCR duplicates, so read 2 shouldn't have a different sequence (unless there are sequencing errors).

If you look up the library prep protocol, you can see that cDNA amplification happens before fragmentation.  Therefore, you can have multiple reads starting at different fragmentation points which nevertheless originate from the same mRNA molecule.  It's also quite plain from looking at the UMI counts output and bam output that all reads with the same gene ID, UMI and cell barcode are marked as duplicates save one, even if the reads aligned to different places in the transcript, and that only the one not marked as a duplicate is being counted in the final UMI count

Thanks. I Just had another look and this appears to be indeed the case, also for other drop-seq protocols apart from 10x. To however accept this as an answer, the only thing missing is if DropSeq-tools is doing the right thing by default or if this can be adjusted that way? The docu is surprisingly unclear about this important point.

It's been a while since I looked at DropSeq data, or its library prep protocol, but one has to presume that the Drop-Seq cookbook protocol handles DropSeq data correctly.  The description cited does seem to indicate that all the reads aligning to a gene are assessed for UMIs, not just all the reads aligning to a particular position.  However, I have in the past run the same Drop-Seq data set through the cookbook and (with some manipulation) through 10XGenomics cellranger, and I got very nearly the same results.

To confirm if this is correct, I've just submitted an issue

This is likely correct also for DropSeq tools and DropSeqPipe but this source is specifically for 10x Genomic Chromium platform. The source code for UMI processing for DropSeq tools can be found here: https://github.com/broadinstitute/Drop-seq

The code does EXACTLY what it says in the first quote.
Let's say you have 5 reads on a cell barcode that align to gene A.  Those reads have the following molecular barcodes:
AAAA,BBBB,AAAA,CCCC,DDDD.
There are 4 UNIQUE molecular barcodes.
The only complication to this is that we collapse barcodes with an edit distance of 1, which takes care of over-inflation of UMI counts due to PCR/Sequencing errors.  If we added a barcode with the sequence AAAB to the original list, it would be collapsed into AAAA as it has an edit distance of 1, and you'd still have 4 unique molecular barcodes.

In retrospect, the quote makes now sense, I was a bit confused what assembled exactly meant in this context, also the additional rational from swbarnes2 about the fragmentation taken place after PCR now explains why it is done this way (hence the accepted answer). Maybe this could be included in the documentation as explanation?


text:
Convert BAM to properly paired FASTQ files
I thought I had figured this one out. But apparently not.
I need to convert a BAM file of paired-end alignments to two FASTQ files of paired reads to realign them, with a twist: I only want reads that fall within a defined region.
Currently I am using the following command (the region is an example; it also fails on different regions; at any rate, this is on hs37d5 but presumably the issue would exist for other references too):
samtools view -b -u -@ 48 alignments.bam 21:31,831,502-31,896,094 \
| samtools collate -n 128 -u -O -@ 48 - /tmp/bam-to-fastq- \
| samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

(Apart from the region filtering, this command is essentially the same that word of god tells us to use.)
The resulting files are not properly paired:
⟩⟩⟩ diff <(gunzip -c reads_R1.fastq.gz | sed -n -e 1~4p -e 10q) \
1 ⟩   <(gunzip -c reads_R2.fastq.gz | sed -n -e 1~4p -e 10q)
1,2d0
< @A00346:17:H52H7DSXX:3:2142:32316:24518
< @A00346:17:H52H7DSXX:3:1153:6705:26412
3a2,3
> @A00346:17:H52H7DSXX:3:1173:9118:3129
> @A00346:17:H52H7DSXX:3:1109:29414:11757

As far as I can see from performing a spot check, these are all reads which are marked as paired in the BAM file, but for which only one mate falls within the filtered region.
And running bwa mem on the files gives me the error:

[mem_sam_pe] paired reads have different names: […]

How can I ensure that only properly paired reads are included in the FASTQ files?

I'm guessing this is just that some reads fall in the target region but their pairs fall outside it. I don't know how to do this (hence a comment), but have you tried extracting only those reads that fall entirely within the target region?

@terdon Right, that’s precisely the problem. However, I’d naively expect samtools fastq to do the right thing here, and actually output pairs of reads, or no read at all. Because resolving this efficiently separately is pretty hard.

The samtools devs might not agree that this is a bug, but the use case you describe is certainly sufficiently common that they should consider implementing it as an option. There's a strong argument that it should really be the default behavior, but that would probably require a major version bump, which raises another set of practical considerations.

That said, I can see some complications with its implementation. We typically see paired reads mapped within a few hundred base pairs, but how will it affect performance if we have very large insert libraries or pairs that are discordantly mapped? In the absence of a data structure that explicitly tracks the location of a read's pair, there are some edge cases that wouldn't be straightforward to handle. From this perspective, I'm not surprised the use case you describe doesn't seem to be supported out-of-the-box.

@DanielStandage samtools collate should take care of this issue: If both reads of a pair are present, they will be adjacent; so no fancy data structure is necessary. In principle this should even be possible using a simple additional pass using sed but I’m currently failing at that.

@KonradRudolph The range-based query (presumably?) requires position-sorted reads, where pairs are separated. Yes, the collate operation should bring the pairs together, but only if both pairs were in the query range. My point was that when reads are sorted by position it's not always straightforward to extend the search to grab pairs of all reads in the query region.

@KonradRudolph Collating before executing the range query would address the issue, but I don't think it's possible to execute a range query on alignments that aren't sorted by genomic position.

@DanielStandage I don’t understand your point then: Filtering of paired reads only needs to happen in the last step, i.e. samtools fastq (or immediately before). That’s why the “current” order of operations (filtering for region, name collation, filtering by proper pairing TBD) is fine, and there’s no issue requiring fancy data structures.

@KonradRudolph I think I misunderstood. I was talking about the difficulty of implementing a strategy that makes sure that if either end of the read is in the query region, then both ends are included in the output. It's also possible to maintain proper pairing by discarding any read whose pair doesn't fall within the query region. After re-reading your post it looks like this is what you were suggesting. And indeed should be easier to implement.

As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of samtools collate.
… | samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

Your samtools fastq command is not doing anything to siphon off these singleton reads. You should add -s /dev/null to get rid of them.
See the discussion in samtools/samtools#874 (especially the part starting here). I have not looked at how the recent changes from jkbonfield, which will land in an upcoming samtools version (after 1.9), might affect this.

Oooh, that explains why I previously thought this worked! My script used to include -s /dev/null but then I removed it because I thought it had no effect. The documentation of these different outputs is somewhat confusing.


text:
R error: could not find function "pathwayAnalysis"
## ------------------------------------------------------------------------
library(limma)
library(statmod)
library(GEOquery)
#library(GEOmetadb)
require(Biobase)
library(ggplot2)
library(rio)
library(affy)
library(affyPLM)
library(reshape)
library(stringr)
library(cowplot)
#library(makecdfenv)
#library(aroma.affymetrix)
#library(RTN)
library(gProfileR)
library(corpcor)
library(viper)
library(rio)
library(dplyr)
library(igraph)
library(RColorBrewer)
library(reshape2)

## ------------------------------------------------------------------------
gse89988 <- getGEO(filename = "GSE89988_family.soft.gz", GSEMatrix = TRUE)

## ------------------------------------------------------------------------
# define samples
samples_nt <- c('GSM2394292', 'GSM2394303')
samples_pam <- c('GSM2394295', 'GSM2394306')
samples_lps <- c('GSM2394296', 'GSM2394307')

# make table with samples metadata
samples <- c(samples_nt, samples_pam, samples_lps)
metadata <- as.data.frame(samples)
colnames(metadata) <- 'sample_id'

# add group information both with the group acronym as a dummy variable -- which
# will be useful for the DEA; also add color for differentiate groups in plots
metadata$group <- c(rep(0, 2), rep(1, 2), rep(2, 2))
metadata$group_name <- c(rep('NT', 2), rep('PAM', 2), rep('LPS', 2))
metadata$group_color <- c(rep('gray', 2), rep('tomato', 2), rep('dodgerblue', 2))

# get positions in GSE object
match(samples, names(GSMList(gse89988)))

# with the following I confirm that the samples correspond to the expected samples
# and get that expression values are already normalized and log2-transformed
GSMList(gse89988)[[49]]
GSMList(gse89988)[[60]]
GSMList(gse89988)[[52]]
GSMList(gse89988)[[63]]
GSMList(gse89988)[[53]]
GSMList(gse89988)[[64]]

## ------------------------------------------------------------------------
#First, we need to make sure that all of the GSMs are from the same platform
gsmplatforms <- lapply(GSMList(gse89988),function(x) {Meta(x)$platform_id})
head(gsmplatforms)
# from the output of the command above we see that indeed all GSMs derive from
# the same platform "GPL17543"

# We can filter the original GSMList to include only those GSMs with the GPL96 platform and use this list for further processing
gsmlist = Filter(function(gsm) {Meta(gsm)$platform_id=='GPL17543'},GSMList(gse89988))
length(gsmlist)

# get the probeset ordering
probesets <- Table(GPLList(gse89988)[[1]])$ID

# make the data matrix from the VALUE columns from each GSM
# being careful to match the order of the probesets in the platform
# with those in the GSMs
data.matrix <- do.call('cbind', lapply(gsmlist,function(x)
                                      {tab <- Table(x)
                                      mymatch <- match(probesets,tab$ID_REF)
                                      return(tab$VALUE[mymatch])
                                     }))
data.matrix <- apply(data.matrix, 2, function(x) {as.numeric(as.character(x))})

# go through the necessary steps to make a compliant ExpressionSet
# and filter for the samples to be analyzed
rownames(data.matrix) <- probesets
colnames(data.matrix) <- names(gsmlist)
data.matrix <- data.matrix[, samples]
pdata <- data.frame(samples = samples)
rownames(pdata) <- samples
pheno <- as(pdata,"AnnotatedDataFrame")
eset <- new('ExpressionSet', exprs = data.matrix, phenoData = pheno)

## ------------------------------------------------------------------------
load_gene_annotation <- function() {

  ifile <- 'GPL17543-38936.txt'
  df <- read.delim(ifile, comment.char = '#')

  return(df)

}

gene_annotation <- load_gene_annotation()

## ------------------------------------------------------------------------
apply_dea <- function(eset, fc_cutoff = 2, adj_p_cutoff = 0.05) {

  # output directory
  ODIR <- 'results/step1/'
  dir.create(ODIR, showWarnings = FALSE, recursive = TRUE)

  # make design matrix
  design <- model.matrix(~ 0 + factor(c(0, 0, 1, 1, 2, 2)))
  colnames(design) <- c("group0", "group1", "group2")

  # plot with the normalized array weights
  aWeights <- arrayWeights(eset, design = design)
  ofile <- paste0(ODIR, 'array_weights_gse89988_samples_used_in_dea.pdf')
  pdf(ofile, width = 10, height = 5)
  dotchart(aWeights, labels = samples,
           groups = as.factor(metadata$group_name),
           pch = 19,
           color = metadata$group_color,
           main = "Array Weights")
  dev.off()

  # linear model fit
  fit <- lmFit(eset, design)
  contrast.matrix <- makeContrasts(group1-group0, group2-group0,
                                   levels = design)
  fit2 <- contrasts.fit(fit, contrast.matrix)
  fit2 <- eBayes(fit2)

  # save
  ODIR='deliverables/'
  dir.create(ODIR, showWarnings = FALSE)
  ofile=paste0(ODIR, 'fit_for_figures_7_and_8.rds')
  saveRDS(fit2, ofile)

  # extract results, annotate and visualize DEA results

  extract_annotate_and_viz_dea_results <- function(coef_index, ref, treat) {

    # Extract DEA results for a given comparison
    df <- topTable(fit2, coef = coef_index, adjust = "BH", number = Inf)
    df$reference <- ref
    df$treatment <- treat

    # remove genes with missing values
    df <- df[complete.cases(df), ]

    # annotate with gene information -- gene annotation has many fields which
    # are not relevant here so I will make a reduced gene information version
    fields_included <- c('ID', 'Species', 'Source', 'Entrez_Gene_ID', 'Symbol')
    df <- merge(df, gene_annotation[fields_included], by.x = 0, by.y = 'ID',
                sort = FALSE)
    colnames(df)[1] <- 'probe_id'

    # prepare data for visualization through volcano plot
    cond1 <- df$adj.P.Val < adj_p_cutoff
    cond2d <- df$logFC < -log2(fc_cutoff)
    cond2u <- df$logFC > log2(fc_cutoff)
    log10p <- -log10(df$adj.P.Val)
    down <- nrow(df[cond1 & cond2d, ])
    up <- nrow(df[cond1 & cond2u, ])
    df$category <- ifelse(cond1 & cond2d, 'down-regulated',
                          ifelse(cond1 & cond2u, 'up-regulated', 'non-regulated'))
    abs_max_logfc <- max(abs(df$logFC))

    # volcano plot
    g <- ggplot(df, aes(x = logFC, y = -log10(adj.P.Val), colour = category)) +
          geom_point(alpha = 0.4, size = 1.75) +
          scale_color_manual(values = c('dodgerblue', 'gray', 'tomato')) +
          xlim(-abs_max_logfc, abs_max_logfc) +
          ylim(c(0, max(-log10(df$adj.P.Val)))) +
          xlab("log2 fold change") +
          ylab("-log10(adjusted p-value)") +
          annotate("text", x = -4, y = 3, size =3,
                   label = paste0(down, "\ndown-regulated DEGs"),
                   color = 'dodgerblue') +
          annotate("text", x = 4, y = 3, size = 3,
                   label = paste0(up, "\nup-regulated DEGs"),
                   color = 'tomato') +
          annotate("text", x = -abs_max_logfc, y = 0, size = 2,
                   label = paste0("DEGs defined with:\nAdj. p-value < ",
                                  adj_p_cutoff,
                                  '\nFold change > ',
                                  fc_cutoff),
                   color = 'black', hjust = 0) +
          ggtitle(paste0(ref, ' vs. ', treat)) +
          theme(plot.title = element_text(hjust = 0.5))

    # save plot to file
    ofile <- paste0(ODIR, 'volcano_plot_', tolower(ref), '_vs_', tolower(treat),
                    '.pdf')
    ggsave(ofile, width = 7, height = 5)

    return(df)

  }

  res <- rbind(extract_annotate_and_viz_dea_results(1, 'NT', 'PAM'),
               extract_annotate_and_viz_dea_results(2, 'NT', 'LPS'))

  # export table with results
  ofile <- paste0(ODIR, 'differential_expression_analysis_results_combined.xlsx')
  export(res, ofile)

}

## ------------------------------------------------------------------------
apply_dea(eset)

After all these i have to run this
gsca <- pathwayAnalysis(fit2, pvaluecutoff=0.05),where fit2 is made after running apply_dea(eset)
As I run the above code i get this

Error in pathwayAnalysis(fit2, pvaluecutoff = 0.05) :    could not
  find function "pathwayAnalysis"

I would share the .rdata file for reproduciblty
rdata
The gse & the annotation files
I would be really glad if i get to know what am i doing wrong.

What package is pathwayAnalysis a part of? Or is it a custom function you wrote?

Note: Sharing RDA file and custom annotation files does not make our jobs easier. It's not fair to expect us to download and use large datasets just so we can solve your problem.

Please, before posing another insanely long question, read this page. The reproducible example should be minimal.

that is absolutely correct regarding large files , regarding  pathwayAnalysis i m actually trying to replicate  this study https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1002223 the codes is part of its supplementary data.i have loaded all the library mentioned but im not sure what is going wrong

@krushnachChandra ha, then my best guess is that they forgot to mention this library.

I made it an answer.

It seems like you are missing function "pathwayAnalysis" presumably from package BLMA, i.e. library('BLMA') should resolve your error once you install that package.
-- edit 1 --
Ha, stupid me. Package BLMA contains file "pathwayAnalysis.R", not the function. A package that contains a FUNCTION called pathwayAnalysis is IntClust. You should try that one.
-- edit 2 --
Alright, stupid me again, now I have not checked that IntClust has capitalized PathwayAnalysis, you need pathwayAnalysis. I found one more function matching the name form package metaX. However, it could be that it's function defined somewhere in their code. You can try to search the source code for example using grep or you try to use google to find packages that contain the function of this name (this is exactly what I was doing, many you should be a bit more careful than I was).

let me try and see if its working

even after installing and loading the library it doesn;t work may it there is something else. I will look into it

Do you get a new error or the same one? (if it's a new one, edit the question and add there the new error to the end)

the same exact error . "Error in pathwayAnalysis(fit2, pvaluecutoff = 0.05) :
could not find function "pathwayAnalysis""

My bad, I was googling too quickly. Check the edit to my answer.

let me give this a try will update ,

and it still doesn;t work same error "  could not find function "pathwayAnalysis"
" after loading that library .May be will have to bit more into the supplementary codes.

Check the new edit. I made another mistake in my quick search, so I have decided to rather give you hints how to resolve your problems by your own (grepping and googling). Hope it help, my time runs out.

will update you

couldn't make it work but anyways thank you for the help ,meanwhile i have written to the paper authors


text:
Subset a protein fasta based on sequence length?
I have a strange protein fasta file. There are several entries for the same gene and I need to extract the entry with the longest width.
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 a
lpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..8944
8))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFGSTGFIGRYVANKLGKIGTQL
ILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKALKYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAG
VEKFIHLSALNCDGPNEAIFSRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWGYKRYDMKYDMFFKLKVSLTQ
KFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLRRMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017771299.1_10 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017771299.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017772069.1_11 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X2] [protein_id=XP_017772069.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109546..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SVTENYFGRIVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQKILKKR
KCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPC
SLLACPE
>lcl|NW_017095468.1_prot_XP_017772815.1_12 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X3] [protein_id=XP_017772815.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD

I've tried using BioStrings but I'm unfamiliar with how to index this type of file. Any help appreciated.

Can you give us an example that includes a duplicated gene? How is the "gene" defined? Should we look for the RefSeq protein ID (XP_123), for the gene=LOC string? All of the above? Will the entire header line be identical or only parts of it?

I've added more of the file to my question. I think the unique/isoform bits are the XP_ (which has information about isoforms) ... NW_ bits are the genes.

NW_017095468 is a whole genome shotgun sequence, it will have many genes. The XP bits are RefSeq protein accessions. These can be different for different protein products of the same gene. It looks like the solutions will need to look at either the gene= or the GeneID: values.

Perhaps more important is for R-MASHup to explain what is "strange" about these sequences and why. It is well known that a gene may encode multiple proteins, so the fact that there are multiple entries per gene doesn't surprise me.

Yes. And choosing the longest is kind of arbitrary, but all this depends on what you're doing. For instance, I once needed to do something very similar because I was looking at various features of certain genes and I wanted to see if protein product length correlated with some other characteristic, so I chose the longest isoform of each gene's products. But, @R-MASHup, if you think the longest will be somehow special (most highly expressed or more representative or anything like that), you will be surprised.

Well it's a nonmodel organism, trying to find orthologs, and so need to reduce the number of entries somehow (esp as it's computationally so slow) ... how would you pick?

For orthologs, I'd probably not pick since it's impossible to predict which protein will have ended up shared with the other organism. I'd just run a tblastn or better exonerate on using all proteins as queries. But yeah, picking the longest might make sense if you don't have a decent machine. Also, if my answer isn't working for you, you probably have ID lines that are different from what you've shown here. My solution just takes the 2nd field of the line ([gene=LOCNNNN] in your example), if that changes, it will not work.

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do:
FastaToTbl file.fa |
    awk '{
            len=length($NF);
            if(len > l[$2]){
                l[$2]=len;
                seq[$2]=$0;
            }
          }
          END{
            for(s in seq){
                print seq[s]
            }
        }' | TblToFasta

Running this on your example gives:
$ FastaToTbl file.fa | awk '{len=length($NF); if(len > l[$2]){l[$2]=len; seq[$2]=$0; }} END{for(s in seq){print seq[s]}}' | TblToFasta
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNI
FLRSNGEYQAKNKFKDFLLLVEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVM
QTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTDDGPLSKILCEEVQENEDTQP
NMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMS
LVSMIDPSPDWIVGVSGLELCLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHP
IRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEKNCVAQVDVSEEDGGVEGDKC
EMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERE
CGPSDRRPLQESKECEGRAWSNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDT
EDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDNSAVVCGFPRDEGGCMSNVDR
WYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGF
STKHRFIRVHPSNGGKRCPQKILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCV
MTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..89448))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFG
STGFIGRYVANKLGKIGTQLILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKAL
KYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAGVEKFIHLSALNCDGPNEAIF
SRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWG
YKRYDMKYDMFFKLKVSLTQKFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLR
RMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS


text:
How to output only the overlapping region of two files?
I have two files of genomic coordinates. I am trying to output only the overlapping region of fileA if it intersects with fileB. Not the original fileA coordinates.
fileA
chr1 25 50
chr1 75 200

fileB
chr1 10 60
chr1 80 90

desired output
chr1 25 50
chr1 80 90

Please provide a small reproducible example using dput so we can help you. All you need to do is run dput(fileA) in R and, paste the output in the terminal.

bedtools is my go-to program for operations on genomic intervals. In particular, the bedtools intersect operation is what you're looking for here.
$ cat fileA
chr1    25      50
chr1    75      200
$ cat fileB
chr1    10      60
chr1    80      90
$ bedtools intersect -a fileA -b fileB
chr1    25      50
chr1    80      90

Make sure that the values are separated by tabs rather than spaces, or bedtools will have trouble reading the files.


text:
How to combine multiple files into one file?
I have multiple files (n=86000) with one column each and I want to combine them all into one file with 86000 columns.
I tried the following command
Paste *.txt > mean.txt

And as a for loop
for i in *.txt, do paste $i > mean.txt; done

This did not appear to work because there are 86000 to combine, I then divided my files into multiple folders with each folder not containing more than 15000 files but this did not working
Format of my single file is as follows (It has 16000 lines)
0.047202943506056
0.0383549703366668
0.0159162908608477
0.0238048058550539
0.0694728928685859
0.0389052658410832

What is the bug in the code?

Thanks. Do you want one long list? In this case the command is 'cat' not 'paste'. Just check you are not  exceeding your disk space.

What is the error message you get?

Equivalently answered here assuming OP wants to concatenate vertically, not horizontally, though the approach could be modified anyway: https://stackoverflow.com/questions/18695105/how-to-concatenate-huge-number-of-files

Re-reading the question: 16000 lines x 86000 files ... this is not alot. I think the error report is important, or else perhaps it just hangs?

@ Fabio Marroni There is no error. When I check the output file it has only one column.

@Michael: No don't want one long list. I want an output file with 86000 columns.

I'm voting to close this question as off-topic because it is not about bioinformatics but file manipulation. In any case it has been answered here.

@Priya The fact that you wanted a file with 86000 columns is not clear from your original post. You simply said you wanted to combine them into a single file. When posting questions, please take time to clearly and precisely describe the issue.

@DanielStandage Apologies for not being clear but I have not got my answer yet.

@Priya Please explain what is not clear from the thread I linked in my first comment.

@DanielStandage Sorry, I didn't click the link in your post. Thank you

@DanielStandage such questions are on topic since manipulating and parsing text files is a very large part of a bioinformatician's work. This question is really not clear, mind you, but it is on topic.

@Priya  tried to clarify the question a bit and I also suspected a typo in the code you tried. I hope the edits I have done are appropriate, if not please roll it back or [edit] it by your own :-).

The bug in your loop appears be "> filename" per column, this over-rides the last iteration resulting in a file with a single column, which has been over-written 85999 times. Thus the file output you observe should be the single column present in the last file of your array. If you replaced "> filename" with ">> filename" within the loop you would then append to the file a new column with each iteration. If the loop fell over at some point you could then simply use one of the outputs of wc filename to assess how many columns you had (you might have to transponse the data for wc to work (lots of scripts for this in Perl), although there will almost certainly be an easy solution for counting columns in awk.
Personally, I would use a dataframe for this operation (possibly R, definately Python pandas), because 86000 columns is a lot of columns in any other walk of life, even if it is just post-processing. Generally, if you just want bash, awk is a better solution in this case (I don't code with it, but some here certainly do). However, I do get the idea that paste will work if each column is an identical length.
BTW paste must be lower case, you've got Paste in the first line of your code, but its okay in the second bit of code.
BTW BTW, I did not know shell would take *.txt directly as an array, with Perl/Python you need to "glob" it and assign it as an array/list.


text:
How to run 'join' command for multiple files (of 2 types) in a folder
I have two types of files for a gene in a folder
File 1: FOS.tf.txt
ABL2
ACTN4
ADGRE5
ADIPOR2
ADRB2
ERCC4
EZR
FAS
FMO4

File 2: FOS.tt.txt
ABL2
CASP9
CCK
DDIT3
ERCC4
EZR
FAS
FMO4

After running command
join FOS.tf.txt FOS.tt.txt > FOS_final.txt

I got output file: FOS_final.txt
ABL2
ERCC4
EZR
FAS
FMO4

My question is that I have these two type of files (tf and tt) for 605 genes in a folder and I want to run join command together for all my genes and get 605 output files.

If all of your .tf.txt and .tt.txt files are in the same directory, the following command should work.
for tffile in $(ls *.tf.txt)
do
    prefix=$(basename $tffile .tf.txt)
    join ${prefix}.tf.txt ${prefix}.tt.txt > ${prefix}_final.txt
done

The basename call can be avoided by anchoring a regex with parameter expansion like so: join ${tffile} ${tffile/%.tf.txt/.tt.txt} > ${tffile/%.tf.txt/_final.txt}

@RamRS Cool! I guess it's up to the individual to decide which ugly solution is less ugly! :-)

Indeed. I wonder if $() is more performant than ${//}.

@ Ram RS:
I am getting this error
${tffile/%.tf.txt/_final.txt}: ambiguous redirect


text:
How can I transform a mapped BAM file into an unmapped BAM file?
In order to use MergeBamAlignment (Picard), I need an unmapped BAM file and a mapped BAM file. I have two mapped BAM files:

one with reads mapped to the reference I want but without metadata such as read group
a second one with reads mapped to the incorrect reference but with metadata

When I try to use MergeBamAlignment with these two files, I get errors, such as "Second read from pair not found in unmapped bam". So I am trying to unmap the reads from my second BAM file. How can I do that? Is there a tool for it? Otherwise, I'll have to manually change the bitwise FLAGs, reference names, CIGARS, etc.

Looks like you need Picard's RevertSam operation (described in context here, manual here).
java -jar picard.jar RevertSam \
     I=input.bam \
     O=reverted.bam

What happens if you use the bam with the metadata as the input to your aligner?  Maybe it will keep all those fields.

bwa mem requires a fastq file, it does not accept a BAM

Okay, I didn't know you had to use that.  You see that bwa aln will take a bam?

good point, the RevertSam worked though


text:
Cluster is split in 2-3 locations on tsne plot - Suerat
I am running a single cell dataset (count data - exon) through Seurat. After running tsne I see a cluster (13) split in 3 different locations on the plot. Here are the commands I am running:
library('Seurat')
## rows (genes) that have 99% zeros were filtered out,
## and then columns (cells) that have 99% zeros were filtered out
isp.data <- read.table(file = "genes.cells99.percent.filtered", sep =
            ",", row.names = 1, header=TRUE)
isp <- CreateSeuratObject(raw.data = isp.data, project =
       "brain_regions")
isp <- NormalizeData(object = isp, normalization.method =
       "LogNormalize",scale.factor = 10000)
isp <- FindVariableGenes(object = isp, mean.function = ExpMean,
       dispersion.function = LogVMR,do.plot = FALSE)
hv.genes <- head(rownames(isp@hvg.info), 1000)
isp <- ScaleData(object = isp, genes.use = hv.genes, num.cores = 3)
isp <- RunPCA(object = isp, pc.genes = hv.genes, pcs.compute = 100,
       do.print = TRUE, pcs.print = 1:5, genes.print = 5)
PCElbowPlot(object = isp, num.pc = 100)
isp <- JackStraw(object = isp, num.replicate = 100, display.progress =
                 FALSE, num.pc = 30)
JackStrawPlot(object = isp, PCs = 1:30)
isp <- FindClusters(object = isp, reduction.type = "pca", dims.use =
        1:26, resolution= 0.75, print.output = 0, save.SNN = TRUE,
        force.recalc=TRUE)
isp <- RunTSNE(isp, dims.use = 1:26, do.fast = T)
TSNEPlot(isp, do.label = T, pt.size = 0.5)

Here is the tsne plot:

I have tried varying parameters for resolution in FindClusters (from 0.5-3.5), but that just increases the clusters found. Further I tried to play around with parameters min.dist (0.75,1), n_neighbors (10 to 30) and perplexity(5 to 50) but I still get at least 1-2 clusters split in the tsne plot.
Another thing I have tried is to reduce the dimension (dims.use=10) on which tsne is ran but the clusters seem to be too close together, like so:

Note: Dimensionality of dataset (after filtering) - 23824 genes across 8967 samples.
My question being why does such behavior happen? Is my filtering strategy not correct? Have I not normalized/scaled data correctly?

Your cluster labels come from graph clustering implemented in the FindClusters() function. The resulting clusters are then visualised with a 2D tSNE plot (via RunTSNE() and TSNEPlot()). So, your cluster 13 is not split into three sub-clusters but cells within cluster 13 look somewhat distant from each other on the tSNE plot. Having mentioned distances on a tSNE plot, please see this SE Data Science entry Can closer points be considered more similar in T-SNE visualization? and especially the comment "What I mean is that you can't just use distance in the lower space as a similarity criterion". In your case "the lower space" would be your tSNE plot.
As far as I can see from your code, you are following the widely-adopted "Seurat workflow" and all seem fine including your TSNE plot output. I could not spot the QC step though: apart from the one at object initiation, one usually applies filters for number of genes and % mitochondrial gene contribution (using FilterCells() in Seurat v2 and subset() in Seurat v3).
It would be a good idea to try UMAP instead of tSNE, please see Evaluation of UMAP as an alternative to t-SNE for single-cell data (Becht et al, 2019). Seurat v3 vignettes are now using UMAP and not tSNE any more.
Lastly, please see How to Use t-SNE Effectively to learn about points to take into account when using/interpreting tSNE plots. You will see that you would get surprisingly dissimilar tSNE plots from the same underlying data with different parameters (some are hardcoded in the RunTSNE() function).

Yes, I was missing step where I needed to filter for mitochondria and number of genes.

Great and very comprehensive answer!

Hi user2998764, I was just curious, do the tSNE plot looks "better" after the QC for % mito and genes?