轉(zhuǎn)錄組分析實(shí)戰(zhàn)附錄：Trinity 拼接結(jié)果質(zhì)量控制

在第二節(jié)中，我們采用了Trinity工具做了轉(zhuǎn)錄組數(shù)據(jù)的拼接分飞，我一共是6個(gè)樣本6個(gè)G的數(shù)據(jù)量肮雨，在我那個(gè)設(shè)置下跑了接近30多個(gè)小時(shí)就完成了拼接工作蕊唐。

那么今天的工作就是通過RSeQC這個(gè)軟件對(duì)拼接結(jié)果進(jìn)行一個(gè)質(zhì)量控制與可視化

這個(gè)軟件主要是針對(duì)于一些臨床RNAseq的數(shù)據(jù)以及有參考基因組的數(shù)據(jù)，但是對(duì)沒有參考基因組的RNAseq數(shù)據(jù)就很多Tool沒有辦法使用简十。

首先，通過bowtie2對(duì)得到的Trinity拼接好的fasta格式進(jìn)行構(gòu)建Index

yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ ll | sort -nk 7total 86667684
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 left.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 right.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:08 both.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:19 .jellyfish_count.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:24 .jellyfish_dump.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:26 .jellyfish_histo.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm_renamed.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 inchworm.K25.L25.DS.fa.finished
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 partitioned_reads.files.list.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 recursive_trinity.cmds.ok
-rw-rw-r-- 1 yeyt yeyt           9 Sep 14 01:08 both.fa.read_count
-rw-rw-r-- 1 yeyt yeyt          10 Sep 14 01:59 inchworm.kmer_count
-rw-rw-r-- 1 yeyt yeyt        2757 Sep 14 10:31 pipeliner.3855.cmds
-rw-rw-r-- 1 yeyt yeyt       22843 Sep 14 01:26 jellyfish.kmers.fa.histo
-rw-rw-r-- 1 yeyt yeyt    13366802 Sep 14 10:39 partitioned_reads.files.list
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 14 10:39 recursive_trinity.cmds
-rw-rw-r-- 1 yeyt yeyt   493022300 Sep 14 03:27 inchworm.K25.L25.DS.fa
-rw-rw-r-- 1 yeyt yeyt 11602365066 Sep 14 01:08 both.fa
-rw-rw-r-- 1 yeyt yeyt 26501596675 Sep 14 01:24 jellyfish.kmers.fa
-rw-rw-r-- 1 yeyt yeyt 49568938526 Sep 14 06:38 scaffolding_entries.sam
drwxrwxr-x 2 yeyt yeyt        4096 Sep 14 10:33 chrysalis/
drwxrwxr-x 3 yeyt yeyt        4096 Sep 14 01:03 insilico_read_normalization/
drwxrwxr-x 4 yeyt yeyt        4096 Sep 14 10:39 read_partitions/
-rw-rw-r-- 1 yeyt yeyt           0 Sep 15 20:51 align_stats.txt
-rw-rw-r-- 1 yeyt yeyt          62 Sep 15 20:51 bowtie2.bam
-rw-rw-r-- 1 yeyt yeyt         651 Sep 15 07:03 Trinity.timing
-rw-rw-r-- 1 yeyt yeyt    10213332 Sep 15 07:03 Trinity.fasta.gene_trans_map
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 15 07:03 recursive_trinity.cmds.completed
-rw-rw-r-- 1 yeyt yeyt   244740565 Sep 15 07:03 Trinity.fasta
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2-build Trinity.fasta Trinity.fasta
Settings:
  Output files: "Trinity.fasta.*.bt2"
  Line rate: 6 (line is 64 bytes)
  Lines per side: 1 (side is 64 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Max bucket size: default
  Max bucket size, sqrt multiplier: default
  Max bucket size, len divisor: 4
  Difference-cover sample period: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  Trinity.fasta
Building a SMALL index
Reading reference sizes
  Time reading reference sizes: 00:00:03
Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
...
Exiting Ebwt::buildToDisk()
Returning from initFromVector
Wrote 103828770 bytes to primary EBWT file: Trinity.fasta.rev.1.bt2
Wrote 55572488 bytes to secondary EBWT file: Trinity.fasta.rev.2.bt2
Re-opening _in1 and _in2 as input streams
Returning from Ebwt constructor
Headers:
    len: 222289920
    bwtLen: 222289921
    sz: 55572480
    bwtSz: 55572481
    lineRate: 6
    offRate: 4
    offMask: 0xfffffff0
    ftabChars: 10
    eftabLen: 20
    eftabSz: 80
    ftabLen: 1048577
    ftabSz: 4194308
    offsLen: 13893121
    offsSz: 55572484
    lineSz: 64
    sideSz: 64
    sideBwtSz: 48
    sideBwtLen: 192
    numSides: 1157761
    numLines: 1157761
    ebwtTotLen: 74096704
    ebwtTotSz: 74096704
    color: 0
    reverse: 1
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ ll | sort -nk 7
total 86822504
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 left.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 right.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:08 both.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:19 .jellyfish_count.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:24 .jellyfish_dump.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:26 .jellyfish_histo.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm_renamed.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 inchworm.K25.L25.DS.fa.finished
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 partitioned_reads.files.list.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 recursive_trinity.cmds.ok
-rw-rw-r-- 1 yeyt yeyt           9 Sep 14 01:08 both.fa.read_count
-rw-rw-r-- 1 yeyt yeyt          10 Sep 14 01:59 inchworm.kmer_count
-rw-rw-r-- 1 yeyt yeyt        2757 Sep 14 10:31 pipeliner.3855.cmds
-rw-rw-r-- 1 yeyt yeyt       22843 Sep 14 01:26 jellyfish.kmers.fa.histo
-rw-rw-r-- 1 yeyt yeyt    13366802 Sep 14 10:39 partitioned_reads.files.list
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 14 10:39 recursive_trinity.cmds
-rw-rw-r-- 1 yeyt yeyt   493022300 Sep 14 03:27 inchworm.K25.L25.DS.fa
-rw-rw-r-- 1 yeyt yeyt 11602365066 Sep 14 01:08 both.fa
-rw-rw-r-- 1 yeyt yeyt 26501596675 Sep 14 01:24 jellyfish.kmers.fa
-rw-rw-r-- 1 yeyt yeyt 49568938526 Sep 14 06:38 scaffolding_entries.sam
drwxrwxr-x 2 yeyt yeyt        4096 Sep 14 10:33 chrysalis/
drwxrwxr-x 3 yeyt yeyt        4096 Sep 14 01:03 insilico_read_normalization/
drwxrwxr-x 4 yeyt yeyt        4096 Sep 14 10:39 read_partitions/
-rw-rw-r-- 1 yeyt yeyt           0 Sep 15 20:51 align_stats.txt
-rw-rw-r-- 1 yeyt yeyt          62 Sep 15 20:51 bowtie2.bam
-rw-rw-r-- 1 yeyt yeyt         651 Sep 15 07:03 Trinity.timing
-rw-rw-r-- 1 yeyt yeyt    10213332 Sep 15 07:03 Trinity.fasta.gene_trans_map
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 15 07:03 recursive_trinity.cmds.completed
-rw-rw-r-- 1 yeyt yeyt   244740565 Sep 15 07:03 Trinity.fasta
drwxrwxr-x 3 yeyt yeyt        4096 Sep 15 18:42 ../
drwxrwxr-x 5 yeyt yeyt        4096 Sep 15 20:51 ./
-rw-rw-r-- 1 yeyt yeyt     1984490 Sep 23 13:50 Trinity.fasta.3.bt2
-rw-rw-r-- 1 yeyt yeyt    55572480 Sep 23 13:50 Trinity.fasta.4.bt2
-rw-rw-r-- 1 yeyt yeyt    55572488 Sep 23 14:03 Trinity.fasta.2.bt2
-rw-rw-r-- 1 yeyt yeyt    55572488 Sep 23 14:16 Trinity.fasta.rev.2.bt2
-rw-rw-r-- 1 yeyt yeyt   103828770 Sep 23 14:03 Trinity.fasta.1.bt2
-rw-rw-r-- 1 yeyt yeyt   103828770 Sep 23 14:16 Trinity.fasta.rev.1.bt2

在最后生成的6個(gè)以bt2結(jié)尾的則是Index文件

接下來進(jìn)行Bowtie2回貼并生成sam文件

yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B251_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B251_2.P.fq.gz -S B251.sam

#最后生成的以下文件log：
#回貼B251的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B251_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B251_2.P.fq.gz -S B251.sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
28213701 reads; of these:
  28213701 (100.00%) were paired; of these:
    3865337 (13.70%) aligned concordantly 0 times
    2140365 (7.59%) aligned concordantly exactly 1 time
    22207999 (78.71%) aligned concordantly >1 times
    ----
    3865337 pairs aligned concordantly 0 times; of these:
      134400 (3.48%) aligned discordantly 1 time
    ----
    3730937 pairs aligned 0 times concordantly or discordantly; of these:
      7461874 mates make up the pairs; of these:
        2553395 (34.22%) aligned 0 times
        273693 (3.67%) aligned exactly 1 time
        4634786 (62.11%) aligned >1 times
95.47% overall alignment rate
#回貼B252的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B252_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/B252_2.P.fq.gz -S B252.sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
24423445 reads; of these:
  24423445 (100.00%) were paired; of these:
    2755943 (11.28%) aligned concordantly 0 times
    2003579 (8.20%) aligned concordantly exactly 1 time
    19663923 (80.51%) aligned concordantly >1 times
    ----
    2755943 pairs aligned concordantly 0 times; of these:
      82738 (3.00%) aligned discordantly 1 time
    ----
    2673205 pairs aligned 0 times concordantly or discordantly; of these:
      5346410 mates make up the pairs; of these:
        1943923 (36.36%) aligned 0 times
        258490 (4.83%) aligned exactly 1 time
        3143997 (58.81%) aligned >1 times
96.02% overall alignment rate
#回貼R251的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/R251_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/R251_2.P.fq.gz -S R251sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
24498964 reads; of these:
  24498964 (100.00%) were paired; of these:
    2605874 (10.64%) aligned concordantly 0 times
    2058157 (8.40%) aligned concordantly exactly 1 time
    19834933 (80.96%) aligned concordantly >1 times
    ----
    2605874 pairs aligned concordantly 0 times; of these:
      68645 (2.63%) aligned discordantly 1 time
    ----
    2537229 pairs aligned 0 times concordantly or discordantly; of these:
      5074458 mates make up the pairs; of these:
        1920173 (37.84%) aligned 0 times
        259673 (5.12%) aligned exactly 1 time
        2894612 (57.04%) aligned >1 times
96.08% overall alignment rate
#回貼R252的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/R252_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/R252_2.P.fq.gz -S R252.sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
23929511 reads; of these:
  23929511 (100.00%) were paired; of these:
    3455581 (14.44%) aligned concordantly 0 times
    1770888 (7.40%) aligned concordantly exactly 1 time
    18703042 (78.16%) aligned concordantly >1 times
    ----
    3455581 pairs aligned concordantly 0 times; of these:
      132348 (3.83%) aligned discordantly 1 time
    ----
    3323233 pairs aligned 0 times concordantly or discordantly; of these:
      6646466 mates make up the pairs; of these:
        2061887 (31.02%) aligned 0 times
        216206 (3.25%) aligned exactly 1 time
        4368373 (65.72%) aligned >1 times
95.69% overall alignment rate
#回貼W251的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/W251_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/W251_2.P.fq.gz -S W251.sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
25553075 reads; of these:
  25553075 (100.00%) were paired; of these:
    3705332 (14.50%) aligned concordantly 0 times
    2003416 (7.84%) aligned concordantly exactly 1 time
    19844327 (77.66%) aligned concordantly >1 times
    ----
    3705332 pairs aligned concordantly 0 times; of these:
      163553 (4.41%) aligned discordantly 1 time
    ----
    3541779 pairs aligned 0 times concordantly or discordantly; of these:
      7083558 mates make up the pairs; of these:
        2021254 (28.53%) aligned 0 times
        226959 (3.20%) aligned exactly 1 time
        4835345 (68.26%) aligned >1 times
96.04% overall alignment rate
#回貼W252的雙端測序結(jié)果
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bowtie2 -x Trinity.fasta -1 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/W252_1.P.fq.gz -2 /home/yeyt/biodata/NH160034/NH160034/cleandata/assembly/W252_2.P.fq.gz -S W252.sam
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_US:en",
        LC_ALL = (unset),
        LC_PAPER = "zh_CN.UTF-8",
        LC_ADDRESS = "zh_CN.UTF-8",
        LC_MONETARY = "zh_CN.UTF-8",
        LC_NUMERIC = "zh_CN.UTF-8",
        LC_TELEPHONE = "zh_CN.UTF-8",
        LC_IDENTIFICATION = "zh_CN.UTF-8",
        LC_MEASUREMENT = "zh_CN.UTF-8",
        LC_TIME = "zh_CN.UTF-8",
        LC_NAME = "zh_CN.UTF-8",
        LANG = "en_US.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
24577100 reads; of these:
  24577100 (100.00%) were paired; of these:
    3173490 (12.91%) aligned concordantly 0 times
    1898984 (7.73%) aligned concordantly exactly 1 time
    19504626 (79.36%) aligned concordantly >1 times
    ----
    3173490 pairs aligned concordantly 0 times; of these:
      112017 (3.53%) aligned discordantly 1 time
    ----
    3061473 pairs aligned 0 times concordantly or discordantly; of these:
      6122946 mates make up the pairs; of these:
        2060673 (33.65%) aligned 0 times
        226885 (3.71%) aligned exactly 1 time
        3835388 (62.64%) aligned >1 times
95.81% overall alignment rate

這個(gè)過程比較消耗時(shí)間撬腾，我們于此同時(shí)做個(gè)簡單質(zhì)量控制報(bào)告

yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ $TRINITY_HOME/util/TrinityStats.pl Trinity.fasta > Trinitystats.log 
#輸出到Trinitystats.log文件
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ cat Trinitystats.log 


################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  110851
Total trinity transcripts:  220498
Percent GC: 42.98

########################################
Stats based on ALL transcript contigs:
########################################

   Contig N10: 4369
   Contig N20: 3291
   Contig N30: 2640
   Contig N40: 2183
   Contig N50: 1802

   Median contig length: 542
   Average contig: 1008.13
   Total assembled bases: 222289920


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

   Contig N10: 3997
   Contig N20: 2867
   Contig N30: 2195
   Contig N40: 1663
   Contig N50: 1157

   Median contig length: 364
   Average contig: 686.86
   Total assembled bases: 76139520

解釋一下上面的結(jié)果螟蝙。

首先做一個(gè)概括拼接得到多少個(gè)基因，得到多少個(gè)轉(zhuǎn)錄本
然后平均的GC含量是多少
接下來做一個(gè)兩個(gè)工作
一個(gè)是基于所有轉(zhuǎn)錄本的contig統(tǒng)計(jì)
一個(gè)是基于所有基因的統(tǒng)計(jì)
N50代表的是

接下來我們將把得到的sam結(jié)果轉(zhuǎn)化成bam結(jié)果并進(jìn)行排序以提供后期的分析文件

yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ ls *sam | grep '25' |xargs -I [] echo 'samtools view -bS [] | samtools sort -o [].sorted.bam ' > samtoolssort.sh
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ cat samtoolssort.sh 
samtools view -bS B251.sam | samtools sort -o B251.sam.sorted.bam 
samtools view -bS B252.sam | samtools sort -o B252.sam.sorted.bam 
samtools view -bS R251sam | samtools sort -o R251sam.sorted.bam 
samtools view -bS R252.sam | samtools sort -o R252.sam.sorted.bam 
samtools view -bS W251.sam | samtools sort -o W251.sam.sorted.bam 
samtools view -bS W252.sam | samtools sort -o W252.sam.sorted.bam
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bash samtoolssort.sh 
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ bash samtoolssort.sh                                    [bam_sort_core] merging from 41 files...
[bam_sort_core] merging from 36 files...
[bam_sort_core] merging from 36 files...
[bam_sort_core] merging from 35 files...
[bam_sort_core] merging from 38 files...
[bam_sort_core] merging from 36 files...
yeyt@ubuntu:~/biodata/NH160034/NH160034/cleandata/assembly/trinity_out_dir$ ll | sort -nk 7
total 234179528
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 left.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:07 right.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:08 both.fa.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:19 .jellyfish_count.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:24 .jellyfish_dump.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 01:26 .jellyfish_histo.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 .iworm_renamed.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 03:32 inchworm.K25.L25.DS.fa.finished
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 partitioned_reads.files.list.ok
-rw-rw-r-- 1 yeyt yeyt           0 Sep 14 10:39 recursive_trinity.cmds.ok
-rw-rw-r-- 1 yeyt yeyt           9 Sep 14 01:08 both.fa.read_count
-rw-rw-r-- 1 yeyt yeyt          10 Sep 14 01:59 inchworm.kmer_count
-rw-rw-r-- 1 yeyt yeyt        2757 Sep 14 10:31 pipeliner.3855.cmds
-rw-rw-r-- 1 yeyt yeyt       22843 Sep 14 01:26 jellyfish.kmers.fa.histo
-rw-rw-r-- 1 yeyt yeyt    13366802 Sep 14 10:39 partitioned_reads.files.list
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 14 10:39 recursive_trinity.cmds
-rw-rw-r-- 1 yeyt yeyt   493022300 Sep 14 03:27 inchworm.K25.L25.DS.fa
-rw-rw-r-- 1 yeyt yeyt 11602365066 Sep 14 01:08 both.fa
-rw-rw-r-- 1 yeyt yeyt 26501596675 Sep 14 01:24 jellyfish.kmers.fa
-rw-rw-r-- 1 yeyt yeyt 49568938526 Sep 14 06:38 scaffolding_entries.sam
drwxrwxr-x 2 yeyt yeyt        4096 Sep 14 10:33 chrysalis/
drwxrwxr-x 3 yeyt yeyt        4096 Sep 14 01:03 insilico_read_normalization/
drwxrwxr-x 4 yeyt yeyt        4096 Sep 14 10:39 read_partitions/
-rw-rw-r-- 1 yeyt yeyt           0 Sep 15 20:51 align_stats.txt
-rw-rw-r-- 1 yeyt yeyt          62 Sep 15 20:51 bowtie2.bam
-rw-rw-r-- 1 yeyt yeyt         651 Sep 15 07:03 Trinity.timing
-rw-rw-r-- 1 yeyt yeyt    10213332 Sep 15 07:03 Trinity.fasta.gene_trans_map
-rw-rw-r-- 1 yeyt yeyt    47753878 Sep 15 07:03 recursive_trinity.cmds.completed
-rw-rw-r-- 1 yeyt yeyt   244740565 Sep 15 07:03 Trinity.fasta
drwxrwxr-x 3 yeyt yeyt        4096 Sep 15 18:42 ../
-rw-rw-r-- 1 yeyt yeyt         821 Sep 23 15:17 Trinitystats.log
-rw-rw-r-- 1 yeyt yeyt     1984490 Sep 23 13:50 Trinity.fasta.3.bt2
-rw-rw-r-- 1 yeyt yeyt    55572480 Sep 23 13:50 Trinity.fasta.4.bt2
-rw-rw-r-- 1 yeyt yeyt    55572488 Sep 23 14:03 Trinity.fasta.2.bt2
-rw-rw-r-- 1 yeyt yeyt    55572488 Sep 23 14:16 Trinity.fasta.rev.2.bt2
-rw-rw-r-- 1 yeyt yeyt   103828770 Sep 23 14:03 Trinity.fasta.1.bt2
-rw-rw-r-- 1 yeyt yeyt   103828770 Sep 23 14:16 Trinity.fasta.rev.1.bt2
-rw-rw-r-- 1 yeyt yeyt         400 Sep 24 13:22 samtoolssort.sh
-rw-rw-r-- 1 yeyt yeyt  3049959975 Sep 24 15:37 R252.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt  3181086895 Sep 24 16:44 W252.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt  3192193677 Sep 24 15:06 R251.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt  3206939510 Sep 24 14:33 B252.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt  3267705730 Sep 24 16:11 W251.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt  3655386513 Sep 24 14:01 B251.sam.sorted.bam
-rw-rw-r-- 1 yeyt yeyt 20770276094 Sep 24 01:49 R252.sam
-rw-rw-r-- 1 yeyt yeyt 21235142607 Sep 24 02:03 B252.sam
-rw-rw-r-- 1 yeyt yeyt 21293400430 Sep 24 02:07 R251sam
-rw-rw-r-- 1 yeyt yeyt 21346715631 Sep 24 02:15 W252.sam
-rw-rw-r-- 1 yeyt yeyt 22197735984 Sep 24 02:29 W251.sam
-rw-rw-r-- 1 yeyt yeyt 24496840308 Sep 24 03:04 B251.sam

這樣我們就得到了6個(gè)sort后的bam文件

采用以下工具

bam_stat.py

clipping_profile.py

inner_distance.py

read_duplication.py

read_GC.py

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末民傻，一起剝皮案震驚了整個(gè)濱河市胰默，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌漓踢，老刑警劉巖牵署，帶你破解...
沈念sama閱讀 217,277評(píng)論 6贊 503
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異喧半，居然都是意外死亡奴迅，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 92,689評(píng)論 3贊 393
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門挺据，熙熙樓的掌柜王于貴愁眉苦臉地迎上來取具，“玉大人，你說我怎么就攤上這事扁耐∠炯欤” “怎么了？”我有些...
開封第一講書人閱讀 163,624評(píng)論 0贊 353
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵婉称，是天一觀的道長块仆。經(jīng)常有香客問我构蹬，道長，這世上最難降的妖魔是什么悔据？我笑而不...
開封第一講書人閱讀 58,356評(píng)論 1贊 293
?港島之戀（遺憾婚禮）
正文為了忘掉前任庄敛，我火速辦了婚禮，結(jié)果婚禮上蜜暑，老公的妹妹穿的比我還像新娘铐姚。我一直安慰自己，他們只是感情好肛捍，可當(dāng)我...
茶點(diǎn)故事閱讀 67,402評(píng)論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布隐绵。她就那樣靜靜地躺著，像睡著了一般拙毫。火紅的嫁衣襯著肌膚如雪依许。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 51,292評(píng)論 1贊 301
城市分裂傳說
那天缀蹄，我揣著相機(jī)與錄音峭跳，去河邊找鬼。笑死缺前，一個(gè)胖子當(dāng)著我的面吹牛蛀醉，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播衅码，決...
沈念sama閱讀 40,135評(píng)論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼拯刁，長吁一口氣：“原來是場噩夢(mèng)啊……” “哼！你這毒婦竟也來了逝段？” 一聲冷哼從身側(cè)響起垛玻，我...
開封第一講書人閱讀 38,992評(píng)論 0贊 275
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎奶躯，沒想到半個(gè)月后帚桩，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 45,429評(píng)論 1贊 314
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡嘹黔，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 37,636評(píng)論 3贊 334
?白月光啟示錄
正文我和宋清朗相戀三年账嚎，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片儡蔓。...
茶點(diǎn)故事閱讀 39,785評(píng)論 1贊 348
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡醉锄，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出浙值，到底是詐尸還是另有隱情恳不，我是刑警寧澤，帶...
沈念sama閱讀 35,492評(píng)論 5贊 345
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布开呐，位于F島的核電站烟勋，受9級(jí)特大地震影響规求，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜卵惦，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,092評(píng)論 3贊 328
男人毒藥：我在死后第九天來索命
文/蒙蒙一阻肿、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧沮尿，春花似錦丛塌、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,723評(píng)論 0贊 22
一樁弒父案赴邻，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至啡捶，卻和暖如春姥敛，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背瞎暑。一陣腳步聲響...
開封第一講書人閱讀 32,858評(píng)論 1贊 269
情欲美人皮
我被黑心中介騙來泰國打工彤敛，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人了赌。一個(gè)月前我還...
沈念sama閱讀 47,891評(píng)論 2贊 370
代替公主和親
正文我出身青樓墨榄，卻偏偏與公主長得像，于是被迫代替她去往敵國和親勿她。傳聞我的和親對(duì)象是個(gè)殘疾皇子袄秩，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 44,713評(píng)論 2贊 354