Creative Commons License
This blog by Tommy Tang is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

My github papge

Friday, July 25, 2014

run a local blast

I was checking the quality of a ChIP-seq fastq with fastqc, and the overrepresented sequence showed that the data contain Illumina Single End Adapter 1 contamination.

I have a file which contains all the possible illumina library sequences in fasta format (got from here), and I manually looked it and found that the single End Adapter 1 sequence did not match the report. I then was wondering  which sequence is the contamination. To do that, I have to do a local blast of  the overrepresented sequence "GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTAGATCGGAAGAGCTCGTATGCC"
on the illumina library sequences.

Luckily, the ncbi_blast was installed on the HPC. One can follow the link here run local blast
first, use the formatdb command to convert the file to a database

$ module load ncbi_blast

# check what arguments are for the command
$ formatdb -

formatdb 2.2.26   arguments:

  -t  Title for database file [String]  Optional
  -i  Input file(s) for formatting [File In]  Optional
  -l  Logfile name: [File Out]  Optional
    default = formatdb.log
  -p  Type of file
         T - protein
         F - nucleotide [T/F]  Optional
    default = T
  -o  Parse options
         T - True: Parse SeqId and create indexes.
         F - False: Do not parse SeqId. Do not create indexes.
 [T/F]  Optional
    default = F
  -a  Input file is database in ASN.1 format (otherwise FASTA is expected)
         T - True,
         F - False.
 [T/F]  Optional
    default = F
  -b  ASN.1 database in binary mode
         T - binary,
         F - text mode.
 [T/F]  Optional
    default = F
  -e  Input is a Seq-entry [T/F]  Optional
    default = F
  -n  Base name for BLAST files [String]  Optional
  -v  Database volume size in millions of letters [Integer]  Optional
    default = 4000
  -s  Create indexes limited only to accessions - sparse [T/F]  Optional
    default = F
  -V  Verbose: check for non-unique string ids in the database [T/F]  Optional
    default = F
  -L  Create an alias file with this name
        use the gifile arg (below) if set to calculate db size
        use the BLAST db specified with -i (above) [File Out]  Optional
  -F  Gifile (file containing list of gi's) [File In]  Optional
  -B  Binary Gifile produced from the Gifile specified above [File Out]  Optional
  -T  Taxid file to set the taxonomy ids in ASN.1 deflines [File In]  Optional


$ formatdb -i adaptors_list.fa -o F -p F

several other files with names ended with nhr, nin and nsq are generated in the same directory.
then run blast:

# check argument for balstall
$ blastall -
blastall 2.2.26   arguments:

  -p  Program Name [String]
  -d  Database [String]
    default = nr
  -i  Query File [File In]
    default = stdin
  -e  Expectation value (E) [Real]
    default = 10.0
  -m  alignment view options:
0 = pairwise,
1 = query-anchored showing identities,
2 = query-anchored no identities,
3 = flat query-anchored, show identities,
4 = flat query-anchored, no identities,
5 = query-anchored no identities and blunt ends,
6 = flat query-anchored, no identities and blunt ends,
7 = XML Blast output,
8 = tabular,
9 tabular with comment lines
10 ASN, text
11 ASN, binary [Integer]
    default = 0
    range from 0 to 11
  -o  BLAST report Output File [File Out]  Optional
    default = stdout
  -F  Filter query sequence (DUST with blastn, SEG with others) [String]
    default = T
  -G  Cost to open a gap (-1 invokes default behavior) [Integer]
    default = -1
  -E  Cost to extend a gap (-1 invokes default behavior) [Integer]
    default = -1
  -X  X dropoff value for gapped alignment (in bits) (zero invokes default behavior)
      blastn 30, megablast 20, tblastx 0, all others 15 [Integer]
    default = 0
  -I  Show GI's in deflines [T/F]
    default = F
  -q  Penalty for a nucleotide mismatch (blastn only) [Integer]
    default = -3
  -r  Reward for a nucleotide match (blastn only) [Integer]
    default = 1
  -v  Number of database sequences to show one-line descriptions for (V) [Integer]
    default = 500
  -b  Number of database sequence to show alignments for (B) [Integer]
    default = 250
  -f  Threshold for extending hits, default if zero
      blastp 11, blastn 0, blastx 12, tblastn 13
      tblastx 13, megablast 0 [Real]
    default = 0
  -g  Perform gapped alignment (not available with tblastx) [T/F]
    default = T
  -Q  Query Genetic code to use [Integer]
    default = 1
  -D  DB Genetic code (for tblast[nx] only) [Integer]
    default = 1
  -a  Number of processors to use [Integer]
    default = 1
  -O  SeqAlign file [File Out]  Optional
  -J  Believe the query defline [T/F]
    default = F
  -M  Matrix [String]
    default = BLOSUM62
..........
..........

$ blastall -i overrepresented_sequence.fa -d adaptors_list.fa -p blastn -e 1e-6 

look at the result, I found that Adapter 2 is the contamination 
>Illumina_Single_End_Apapter_2
          Length = 34

 Score = 63.9 bits (32), Expect = 1e-14
 Identities = 32/32 (100%)
 Strand = Plus / Minus


Query: 1  gatcggaagagctcgtatgccgtcttctgctt 32
          ||||||||||||||||||||||||||||||||
Sbjct: 33 gatcggaagagctcgtatgccgtcttctgctt 2


I then feed Trimmomatic with this fasta file and removed the adaptor.



Tuesday, July 22, 2014

use python to change the header of a fasta file based on a dictionary in another file

I saw on SEQanwser, someone asked this question here, and I used python to resolve this problem.
I have not written in python for so long, I've been using R for a long time. This provide me a good practice opportunity.

The basic problem is to format some files, and I know it can be done by other tools like sed and awk ( see an awk solution below), but python just gives you much more flexibility in terms of what you can do and gives you more readability.
myfasta file

>comp1558_c0_seq1 len=204 path=[4976:0-203]
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 len=222 path=[8835:0-221]
GTACAATCACGGCACATT......
>comp8570_c0_seq1 len=232 path=[3344:0-232]
GCATCGCTTATCGGTTTA......

annotation file:

comp1558_c0_seq1 repressor protein
comp8142_c0_seq1 aspartate aminotransferase
comp8357_c0_seq1 l-xylulose reductase
comp8387_c0_seq1 succinyl- synthetase alpha
comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase

see the code below:




output:
>comp1558_c0_seq1 repressor protein
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 aspartate aminotransferase
GTACAATCACGGCACATT......
>comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase
GCATCGCTTATCGGTTTA......

if the dictionary file has only two records for each line, this awk should work:
http://stackoverflow.com/questions/11678939/replace-text-based-on-a-dictionary
NR == FNR {
  rep[$1] = $2
  next
} 

{
    for (key in rep) {
      gsub(key, rep[key])
    }
    print
}

Friday, July 4, 2014

get an idea of a gene expression value across samples by GEOquery

My next-door PI wanted to look at the TRIM29 expression levels in a series of tumor xenografts from a microarray data. I used GEOquery bioconductor package to get the log2 transformed values and plot a boxplot for him. see the code below. Also, I had a previous post on GEOquery http://crazyhottommy.blogspot.com/2013/12/geoquery-to-access-geo-datasets.html



after got the Expression set object, I looked the names of the object:
> names(pData(gset))
 [1] "title"                    "geo_accession"            "status"                   "submission_date"        
 [5] "last_update_date"         "type"                     "channel_count"            "source_name_ch1"        
 [9] "organism_ch1"             "characteristics_ch1"      "biomaterial_provider_ch1" "molecule_ch1"          
[13] "extract_protocol_ch1"     "label_ch1"                "label_protocol_ch1"       "taxid_ch1"              
[17] "source_name_ch2"          "organism_ch2"             "characteristics_ch2"      "biomaterial_provider_ch2"
[21] "molecule_ch2"             "extract_protocol_ch2"     "label_ch2"                "label_protocol_ch2"    
[25] "taxid_ch2"                "hyb_protocol"             "scan_protocol"            "description"            
[29] "data_processing"          "platform_id"              "contact_name"             "contact_email"          
[33] "contact_phone"            "contact_department"       "contact_institute"        "contact_address"        
[37] "contact_city"             "contact_state"            "contact_zip/postal_code"  "contact_country"        
[41] "supplementary_file"       "data_row_count"

then, only look at the meta data:
> pData(gset)[,c(1,2,28)]
                                 title geo_accession  description
GSM847887                  ML-5998-TG1     GSM847887    Xenograft
GSM847888              Baylor 2147 TG6     GSM847888    Xenograft
GSM847890             Baylor 2665A TG6     GSM847890    Xenograft
GSM847891              Baylor 3107 TG5     GSM847891    Xenograft
GSM847892              Baylor 3143 TG5     GSM847892    Xenograft
GSM847893              Baylor 3204 TG5     GSM847893    Xenograft
GSM847894              Baylor 3561 TG5     GSM847894    Xenograft
GSM847895              Baylor 3611 TG5     GSM847895    Xenograft
GSM847896             Baylor 3613A TG5     GSM847896    Xenograft
GSM847898              Baylor 3807 TG5     GSM847898    Xenograft
GSM847901              Baylor 3887 TG5     GSM847901    Xenograft
GSM847902              Baylor 3904 TG5     GSM847902    Xenograft
GSM847903              Baylor 3936 TG5     GSM847903    Xenograft
GSM847904              Baylor 3963 TG5     GSM847904    Xenograft
GSM847905              Baylor 4013 TG1     GSM847905    Xenograft
GSM847907              Baylor 4175 TG1     GSM847907    Xenograft
GSM847908              Baylor 4195 TG4     GSM847908    Xenograft
GSM847909              Baylor 4272 TG1     GSM847909    Xenograft
GSM847911              Baylor 4664 TG1     GSM847911    Xenograft
GSM847914              Baylor 4888 TG1     GSM847914    Xenograft
GSM847915              Baylor 4913 TG1     GSM847915    Xenograft
GSM847917                  ML-4189-TG2     GSM847917    Xenograft
GSM847919                  ML-5097-TG2     GSM847919    Xenograft
GSM847920                  ML-5156-TG2     GSM847920    Xenograft
GSM847921                  ML-5471-TG2     GSM847921    Xenograft
GSM847922 9830-000060B NEW PROTOCOL V5     GSM847922 Tumor sample
GSM847923            9830-000094B-244K     GSM847923 Tumor sample
GSM847924          9830-000424B-244Kv5     GSM847924 Tumor sample
GSM847925          9830-000517B-244Kv5     GSM847925 Tumor sample
GSM848100 9830-010118B NEW PROTOCOL V5     GSM848100 Tumor sample
GSM848101 9830-010130B NEW PROTOCOL V5     GSM848101 Tumor sample
GSM848102 9830-010214B NEW PROTOCOL V5     GSM848102 Tumor sample
GSM848103 9830-010255B NEW PROTOCOL V5     GSM848103 Tumor sample
GSM848104 9830-010384B NEW PROTOCOL V5     GSM848104 Tumor sample
GSM848105          9830-010461B-244Kv5     GSM848105 Tumor sample
GSM848106        9830-020018B-244K2008     GSM848106 Tumor sample
GSM848107          9830-020025B-244Kv5     GSM848107 Tumor sample
GSM848108          9830-020039B-244Kv5     GSM848108 Tumor sample
GSM848109        9830-020185B-244K2008     GSM848109 Tumor sample
GSM848110          9830-020310B-244Kv5     GSM848110 Tumor sample
GSM848111 9830-020340B NEW PROTOCOL V5     GSM848111 Tumor sample
GSM848112 9830-020416B NEW PROTOCOL V5     GSM848112 Tumor sample
GSM848113        9830-030267B-244K2008     GSM848113 Tumor sample
GSM848114          9830-030446B-244Kv5     GSM848114 Tumor sample
GSM848115        9830-030597B-244K2008     GSM848115 Tumor sample
GSM848116         UNC-000279B-244K2008     GSM848116 Tumor sample
GSM848117         UNC-010208B-244K2008     GSM848117 Tumor sample
GSM848118           UNC-010224B-244Kv5     GSM848118 Tumor sample
GSM848119           UNC-010304B-244Kv5     GSM848119 Tumor sample
GSM848120         UNC-010509B-244K2008     GSM848120 Tumor sample
GSM848121           UNC-020155B-244Kv5     GSM848121 Tumor sample
GSM848122             UNC-020320B-244K     GSM848122 Tumor sample
GSM848123         UNC-020578B-244k-DGE     GSM848123 Tumor sample
GSM848124         UNC-030065B-244K2008     GSM848124 Tumor sample
GSM848125         UNC-030183B-244K2008     GSM848125 Tumor sample
GSM848126         UNC-030370B-244K2008     GSM848126 Tumor sample
GSM848127         UNC-030528B-244K2008     GSM848127 Tumor sample
GSM848128         UNC-040011B-244K2008     GSM848128 Tumor sample
GSM848129           UNC-960028B-244Kv5     GSM848129 Tumor sample
GSM848130           WashU-15720-244Kv5     GSM848130 Tumor sample

there are 10 probes for TRIM29 gene:
> gpl_ann[gpl_ann$"Blast Gene Symbol"=="TRIM29",]
           ID    GB_ACC SPOT_ID Public id       Probe NAME Blast Gene ID Blast Refseq ID Blast Gene Symbol
35972   35972 NM_012101         NM_012101 NM_012101_2_2683         23650     NM_012101.3            TRIM29
39854   39854 NM_012101         NM_012101     A_23_P340123         23650     NM_012101.3            TRIM29
76447   76447 NM_012101         NM_012101     A_23_P340123         23650     NM_012101.3            TRIM29
94361   94361 NM_012101         NM_012101     A_23_P203267         23650     NM_012101.3            TRIM29
97602   97602 NM_012101         NM_012101 NM_012101_2_2608         23650     NM_012101.3            TRIM29
111151 111151 NM_012101         NM_012101 NM_012101_2_2683         23650     NM_012101.3            TRIM29
120477 120477 NM_012101         NM_012101     A_23_P203260         23650     NM_012101.3            TRIM29
132274 132274 NM_012101         NM_012101     A_23_P203260         23650     NM_012101.3            TRIM29
187541 187541 NM_012101         NM_012101 NM_012101_2_2608         23650     NM_012101.3            TRIM29
190531 190531 NM_012101         NM_012101     A_23_P203267         23650     NM_012101.3            TRIM29
               Blast Gene Description Blast Chromosome Map Location
35972  tripartite motif-containing 29                     11q22-q23
39854  tripartite motif-containing 29                     11q22-q23
76447  tripartite motif-containing 29                     11q22-q23
94361  tripartite motif-containing 29                     11q22-q23
97602  tripartite motif-containing 29                     11q22-q23
111151 tripartite motif-containing 29                     11q22-q23
120477 tripartite motif-containing 29                     11q22-q23
132274 tripartite motif-containing 29                     11q22-q23
187541 tripartite motif-containing 29                     11q22-q23
190531 tripartite motif-containing 29                     11q22-q23

the boxplot of the TRIM29 expression levels across different samples: