基因課FTP地址：ftp://http://gsx.genek.tv/2020-3-10%E7%9B%B4%E6%92%AD%E4%B8%80%E4%B8%AA%E5%AE%8C%E6%95%B4%E7%9A%84%E8%BD%AC%E5%BD%95%E7%BB%84%E9%A1%B9%E7%9B%AE/
聽張旭東老師的課

服務(wù)器間數(shù)據(jù)拷貝

兩臺服務(wù)器間的數(shù)據(jù)拷貝 用 scp 用戶名@服務(wù)器：文件路徑

樣本名稱處理

• sed -i 's/.fastq.gz//' xxx.file,類似rename的方式處理文件內(nèi)容
• awk 利用文件名生成樣本信息表

構(gòu)建參考基因組

hisat2-build genome.fasta genome
運行時間 —— 二十分鐘以內(nèi)

比對

hisat2 -x [ref-genome] -U [input_filename].fq.gz -S [output_file name].sam -p [threads] --new-summary --rna-strandness R

• -U 寫基因組名，不寫后綴
• PE 的輸入文件-U項換為 -1 和 -2
• 線程數(shù)建議一般設(shè)置為2-6即可
• --new-summary 歷史原因，使用tophat的舊日志文件格式則不加此項，用新格式日志文件則加此項
• --rna-strandness 鏈特異性文庫
  鏈特異性測序針對性解決的問題是：某些基因所在的正鏈與另一些基因的反鏈有交集，表達量定給誰？
  若不是鏈特異性測序，去掉此項
  若是鏈特異性測序，要問清楚是用的哪個技術(shù)， 大部分都用的是dUTP(90%) 如果是，
  單末端測序(SE) --rna-strandness參數(shù)設(shè)置為R
  PE 設(shè)置為RF
• 目前絕大部分為鏈特異性測序
• 非連特異性測序按照鏈特異性測序比對，有問題
  連特異性測序按照非鏈特異性測序比對，問題不大
• 對于網(wǎng)上下載的測序數(shù)據(jù)有一些沒有寫明是鏈特異性測序還是非鏈特異性測序，解決方法：先假設(shè)是非連特異性測序，比對后在IGV上發(fā)現(xiàn)序列都是同一方向，則為鏈特異性文庫
• 比對DNA序列到基因組用bwa軟件
• 批量生成比對腳本，用awk實現(xiàn)
  vim的Ctrl+v也可以實現(xiàn)
  linux for循環(huán)也可以實現(xiàn)
• 一般不需要設(shè)置錯配率，默認就好。若比對后發(fā)現(xiàn)比對率特別低，則需要考慮。
• 比對率一般至少70以上，比對率和 參考基因組測序組裝質(zhì)量、比對軟件、測序品種與參考基因組物種親緣關(guān)系 相關(guān)
• 并行總線程可超過CPU數(shù)，超過即排隊

比對結(jié)果查看

PE比對結(jié)果log文件

比對率97.44%

比對結(jié)果比對率統(tǒng)計與可視化

比對率結(jié)果在.log文件中
用軟件MultiQC將多個log文件的結(jié)果統(tǒng)計

比對結(jié)果壓縮排序

samtools sort -o xxx.bam xxx.sam

• 這步比較耗內(nèi)存，可以一個一個來
• samtools view是只負責(zé)sam → bam的格式轉(zhuǎn)換
• bam文件構(gòu)建好后，sam文件就可以刪除了

對一個bam文件進行統(tǒng)計

samtools flagstat xxx.bam

• 統(tǒng)計比對率(不同軟件比對出來有差異)
• 不同比對策略算出來的比對率略有不同，有primery的比對和secondary的比對，有區(qū)別
• 轉(zhuǎn)錄組中建議以hisat2統(tǒng)計結(jié)果為準(zhǔn)
  samtools的統(tǒng)計是通用的，沒有對特定軟件進行優(yōu)化

構(gòu)建bam index

samtools index xxx.bam

• 構(gòu)建index在轉(zhuǎn)錄組分析中除了IGV展示，沒有其他用處
• 在重測序中，bam文件構(gòu)建index可用于變異檢測等

IGV可視化

• IGV官網(wǎng) software.broadinstitute.org/software/igv/download
• 用Java寫的軟件，跨平臺(優(yōu)點)，易報錯(劣)
• 需要的文件
  ①導(dǎo)入基因組文件 genome.fasta
  ②基因注釋文件genes.gtf
  ③sample.bam
  ④sample.bam.bai
• 使用步驟
  • 建立基因組庫
    Genomes → Creat .genome File...

  • 加載bam文件
    File → Load from file

    
    
    igv基因組建立
    
    
    gtf文件內(nèi)容差不多長這樣，只有transcript和exon

igv查看mapping結(jié)果

• 桌面軟件
• 同一基因區(qū)同一堿基處，若大概一半一半的A/C概率，則為雜合，若極少量SNP可能是測序錯誤或RNA編輯
• DNA測序鑒定突變：參考基因組為A,測序基因組絕大部分是A
• 個體重測序更關(guān)心的是基因分型，不是只關(guān)心變異區(qū)域，表型=基因+環(huán)境，表型不完全有基因型決定。關(guān)心的是比例問題，概率問題

附加：轉(zhuǎn)錄本基因結(jié)構(gòu)組裝

• hisat2對應(yīng)軟件為stringtie
• 如果別人的基因組組裝做的太差，可選擇自己組裝
• 自己組裝建議使用PASA流程，組裝轉(zhuǎn)錄本

代碼集中營

nohup hisat2-build xxx_genome.fasta xxx_genome 1>hisat2-build.log 2>&1 & # 標(biāo)準(zhǔn)輸出與錯誤輸出到同一文件
# 比對
# SE
hisat2 -x [ref-genome] -U [input_filename].fq.gz -S [output_file name].sam -p [threads] --new-summary --rna-strandness R 1>hisat2.log 2>&1
# PE
hisat2 -x [ref-genome] -1 [input_filename]_1.fq -2 [input_filename]_2.fq -S [output_file name].sam -p [threads] --new-summary --rna-strandness RF 1>hisat2.log 2>&1