劉小澤寫于19.7.30
來北大的龍星計劃做助教，接觸了一些之前沒有了解過的一些知識，簡單做下記錄，我個人其實(shí)是沒有這方面需求的
因?yàn)檎n程都是配置在華為云上的，所以使用的數(shù)據(jù)也都是配置好的，感興趣的小伙伴其實(shí)也不需要自己跑流程(因?yàn)橛玫囊捕际菧y試數(shù)據(jù))，看看怎么使用參數(shù)即可
課程的github 在https://github.com/WGLab/dragonstar2019.git
激活conda
如果看到這里的你有l(wèi)inux基礎(chǔ)，那么就可以自己配置conda，自主下載安裝軟件與環(huán)境了
如果發(fā)現(xiàn)自己賬號下面shared這個目錄沒有內(nèi)容(出現(xiàn)No such file or directory 可以聯(lián)系華為云)
# 教程已經(jīng)將conda封裝好，直接激活使用
. /shared/miniconda3/etc/profile.d/conda.sh
conda activate base
進(jìn)入測試目錄
mkdir -p ~/project/alignment/short-reads
cd ~/project/alignment/short-reads
# 然后軟鏈接數(shù)據(jù)到當(dāng)前目錄(data、ref、vcf)
ln -s /shared/data/NA12878_test_data/short-reads data
ln -s /shared/data/ref_hg37_chr1/ref ref
ln -s /shared/data/ref_hg37_chr1/vcf vcf
以下全部操作都在
~/project/alignment/short-reads目錄下進(jìn)行
質(zhì)控
看下數(shù)據(jù)質(zhì)量
mkdir -p fastqc_output
fastqc -t 2 -o fastqc_output -f fastq data/sr.chr1.2mb_1.fq data/sr.chr1.2mb_2.fq
結(jié)果的Html文件用filezilla導(dǎo)出(選擇自己合適的版本進(jìn)行安裝=》如何使用請看小潔的解答)
(如果不聯(lián)網(wǎng)的話，可能圖片看不了)
對于結(jié)果的解讀：
短序列比對
測試 reads are from a 2MB region in chr1
看下數(shù)據(jù)量
wc -l data/sr.chr1.2mb_1.fq data/sr.chr1.2mb_2.fq
# 2432672 data/sr.chr1.2mb_1.fq
# 2432672 data/sr.chr1.2mb_2.fq
# 4865344 total
# fq文件4行為1個read 因此有1.2 million paired-end reads
方法一 minimap2
這款軟件也是BWA的作者李恒開發(fā)的，它其實(shí)更適合于比對PacBio或Nanopore的DNA/mRNA數(shù)據(jù)
關(guān)于程序怎么用，直接命令行輸入minimap2
minimap2 -ax sr ref/hg37d5.chr1.fa data/sr.chr1.2mb_1.fq data/sr.chr1.2mb_2.fq | samtools sort | samtools view -bS - > chr1.2mb.mp2.bam
# -ax sr 參數(shù)含義： "short genomic paired-end reads"
# 大約1分鐘
# 構(gòu)建索引
samtools index chr1.2mb.mp2.bam
方法二 bwa-mem
bwa mem ref/hg37d5.chr1.fa data/sr.chr1.2mb_1.fq data/sr.chr1.2mb_2.fq | samtools sort | samtools view -bS - > chr1.2mb.bwa.bam
# 大約4min
# 可以加個參數(shù)（-t）加速
bwa mem -t 2 ref/hg37d5.chr1.fa data/sr.chr1.2mb_1.fq data/sr.chr1.2mb_2.fq | samtools sort | samtools view -bS - > chr1.2mb.bwa.bam
# 大約2m3.722s
# 索引
samtools index chr1.2mb.bwa.bam
運(yùn)行命令的同時可以新開一個窗口，登錄進(jìn)去輸入top看看運(yùn)行進(jìn)程，看到用兩個CPU的時候CPU使用量是200%左右
查看bam文件
SAM/BAM文件的格式：https://samtools.github.io/hts-specs/SAMv1.pdf
簡單查看
samtools view chr1.2mb.bwa.bam | less -SN # 按q鍵退出
# 或者直接打印前10行
samtools view chr1.2mb.bwa.bam | head
使用samtools工具
# 格式： samtools tview [options] <aln.bam> [ref.fasta]
# 1.先定義位置，再看
samtools tview -p 1:156000000 chr1.2mb.bwa.bam ref/hg37d5.chr1.fa
# -p chr:pos go directly to this position
# 2.先看，再定義位置
samtools tview chr1.2mb.bwa.bam ref/hg37d5.chr1.fa
# 然后輸入g，就會出現(xiàn)一個對話框，輸入1:156000211跳轉(zhuǎn)(如下圖)
# 退出這個界面，按q
這里就有人問了？為什么我看到的是類似這種
，而且我的還有其他字符，例如既有,，又有.，又有大寫字母，又有小寫字母，應(yīng)該怎么理解，其實(shí)簡單的搜索就可以解決這個問題：我們只需要搜索samtools tview symbols類似的字符串，就能找到答案
例如在：https://www.biostars.org/p/151739/就有解釋：
找變異
主要利用bcftools這個工具，詳見：https://samtools.github.io/bcftools/bcftools.html#mpileup
1 對minimap2產(chǎn)出的bam操作
bcftools mpileup -r 1:155000000-157000000 -f ref/hg37d5.chr1.fa chr1.2mb.mp2.bam | bcftools call -mv -Ob -o mp2.bcftools.call.bcf
# 這一步比較慢，大約需要4m25.991s
關(guān)于bcftools mpileup
- -r: regions
- -f: fasta-ref FILE
關(guān)于bcftools call
-m: multiallelic-caller
-v: variants-only
-
-O: output-type b|u|z|v (詳見：https://samtools.github.io/bcftools/bcftools.html#common_options)
compressed BCF (b); uncompressed BCF (u); compressed VCF (z); uncompressed VCF (v)
那么b和u有什么區(qū)別呢？可以看到上面那個文件是用b參數(shù)，下面那個使用u參數(shù)得到的，大小差別很大
-o：output FILE
然后查看結(jié)果
# 直接看
bcftools view mp2.bcftools.call.bcf | less -SN
# 使用過濾條件(如：只看質(zhì)量值大于20的變異)
bcftools view -i '%QUAL>=20' mp2.bcftools.call.bcf | less -SN
# -i: include EXPRESSION
像這種低質(zhì)量的就會被過濾掉
因?yàn)閎cf文件是二進(jìn)制，有時需要進(jìn)行轉(zhuǎn)換成文本文件vcf
bcftools view mp2.bcftools.call.bcf > mp2.bcftools.call.vcf
2 對bwa mem產(chǎn)出的bam操作
bcftools mpileup -r 1:155000000-157000000 -f ref/hg37d5.chr1.fa chr1.2mb.bwa.bam | bcftools call -mv -Ob -o bwa.bcftools.call.bcf
bcftools view bwa.bcftools.call.bcf > bwa.bcftools.call.vcf
假如說出現(xiàn)這種報錯：
就需要重新運(yùn)行
samtools index chr1.2mb.bwa.bam這個命令
對變異結(jié)果檢驗(yàn)
bcftools view -i '%QUAL>=200' mp2.bcftools.call.bcf > mp2.bcftools.call.vcf
bgzip -f mp2.bcftools.call.vcf -c > mp2.bcftools.call.vcf.gz
bcftools index -f mp2.bcftools.call.vcf.gz
bcftools stats mp2.bcftools.call.vcf.gz | grep "^SN"
# 結(jié)果SN：Summary Numbers
SN 0 number of samples: 1
SN 0 number of records: 2060
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 1765
SN 0 number of MNPs: 0
SN 0 number of indels: 295
SN 0 number of others: 0
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
# 發(fā)現(xiàn)：1765 snps and 295 indels among 2 millions bp region
# 也就是約2000/2M = 1/1000，和理論上人類基因組中有千分之一的突變相符
另外看看transitions/transversions的情況
bcftools stats mp2.bcftools.call.vcf.gz | grep "TSTV"
# TSTV, transitions/transversions:
# TSTV [2]id [3]ts [4]tv [5]ts/tv [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV 0 1260 505 2.50 1260 505 2.50
其中ts是transitions，也就是轉(zhuǎn)換；tv是transversions也就是顛換
轉(zhuǎn)換：嘌呤和嘌呤之間的替換，或嘧啶和嘧啶之間的替換。
顛換：嘌呤和嘧啶之間的替換。
理論上，全基因組分析中Transition/Transversion比值應(yīng)該在2左右，全外顯子分析中比值在3左右
現(xiàn)在將我們得到的和標(biāo)準(zhǔn)的vcf進(jìn)行對比
bcftools isec mp2.bcftools.call.vcf.gz vcf/vcf.chr1.2mb.vcf.gz -p minimap2_perf
# bcftools isec [OPTIONS] A.vcf.gz B.vcf.gz […]
# -p: --prefix DIR
# 結(jié)果會在minimap2_perf目錄中找到四個文件：0000.vcf 0001.vcf 0002.vcf 0003.vcf README.txt
# 其中：0000.vcf 記錄了我們得到的vcf：mp2.bcftools.call.vcf.gz
# 0001.vcf 記錄了標(biāo)準(zhǔn)的vcf：vcf/vcf.chr1.2mb.vcf.gz
# 0002.vcf 記錄了我們得到的vcf中有多少與標(biāo)準(zhǔn)vcf有交集
# 接下來主要用0002.vcf
bcftools stats minimap2_perf/0002.vcf | grep "^SN"
# 結(jié)果如下：
SN 0 number of samples: 1
SN 0 number of records: 1588
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 1569
SN 0 number of MNPs: 0
SN 0 number of indels: 19
SN 0 number of others: 0
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
# 可以看到這里的交集中有1588個變異，也就是我們檢測的變異中有1588個是標(biāo)準(zhǔn)的，之前總共得到了2060個，那么比例就是1588/2060=0.771。如果只考慮SNP的話，就是1569/1765=0.89
關(guān)于
bcftools isec生成的4個文件解釋：最后兩個都是交集為何文件大小差別這么大
https://www.biostars.org/p/224948/ 簡單理解就是：因?yàn)榧词故峭粋€突變，但在不同的vcf包含的注釋信息是不同的
長序列比對
主要使用
minimap2
minimap2 --MD -ax map-ont ref/hg37d5.chr1.fa data/np.chr1.2mb.fq | samtools sort | samtools view -bS - > chr1.2mb.bam
samtools index chr1.2mb.bam
# 關(guān)于x的參數(shù)說明，有以下幾種：
# 1. map-pb/map-ont: PacBio/Nanopore vs reference mapping
# 2. ava-pb/ava-ont: PacBio/Nanopore read overlap
# 3. asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5% sequence divergence
# 4. splice: long-read spliced alignment
# 5. sr: genomic short-read mapping
# 這里的測試數(shù)據(jù)是nanopore的，所以選擇map-ont(ont就是Oxford Nanopore Technologies縮寫)
比對完依然是檢查bam文件
samtools view chr1.2mb.bam | less
samtools tview -p 1:156000000 chr1.2mb.bam ref/hg37d5.chr1.fa
然后找變異
這里不和短序列一樣使用bcftools，而使用了longshot軟件，這是一款2019年放出的軟件，https://github.com/pjedge/longshot 官方介紹其實(shí)很詳細(xì)了
#先激活一個環(huán)境
conda activate longshot
# 它的用法是：longshot [FLAGS] [OPTIONS] --bam <BAM> --ref <FASTA> --out <VCF>
# 然后全部使用默認(rèn)參數(shù)運(yùn)行
longshot --bam chr1.2mb.bam --ref ref/hg37d5.chr1.fa --out mp2.longshot.o.vcf
同樣看看結(jié)果
bcftools view -i '%QUAL>=30' mp2.longshot.o.vcf > mp2.longshot.vcf
bgzip -f mp2.longshot.vcf -c > mp2.longshot.vcf.gz
bcftools index -f mp2.longshot.vcf.gz
bcftools stats mp2.longshot.vcf.gz | grep "TSTV"
# TSTV, transitions/transversions:
# TSTV [2]id [3]ts [4]tv [5]ts/tv [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV 0 1616 529 3.05 1616 529 3.05
看到TS比TV多了3倍
bcftools stats mp2.longshot.vcf.gz | grep "^SN"
# 結(jié)果得到了2145個snp，沒有indel
SN 0 number of samples: 1
SN 0 number of records: 2145
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 2145
SN 0 number of MNPs: 0
SN 0 number of indels: 0
SN 0 number of others: 0
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
與金標(biāo)準(zhǔn)對比
bcftools isec mp2.longshot.vcf.gz vcf/vcf.chr1.2mb.vcf.gz -p longshot_perf
bcftools stats longshot_perf/0002.vcf | grep "^SN"
# 我們得到的結(jié)果與金標(biāo)準(zhǔn)的交叉結(jié)果
SN 0 number of samples: 1
SN 0 number of records: 1189
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 1189
SN 0 number of MNPs: 0
SN 0 number of indels: 0
SN 0 number of others: 0
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
原來短讀長比對金標(biāo)準(zhǔn)結(jié)果總共得到了1588個變異，這里長讀長得到了金標(biāo)準(zhǔn)1189個，計算recall ratio= 1189/1588=0.75；然后我們自己計算得到了2145個，金標(biāo)準(zhǔn)有1189個，計算precision rate = 1189/2145=0.55 ，而這個比例在短讀長里面是0.7以上的，因此長讀長這個方面略顯劣勢
那么短讀長和長讀長之間有多少檢測一致呢？
bcftools isec ../short-reads/mp2.bcftools.call.vcf.gz mp2.longshot.vcf.gz -p shortlong_comparison
bcftools stats shortlong_comparison/0002.vcf | grep "^SN"
# 結(jié)果看到有1294個交集
SN 0 number of samples: 1
SN 0 number of records: 1294
SN 0 number of no-ALTs: 0
SN 0 number of SNPs: 1294
SN 0 number of MNPs: 0
SN 0 number of indels: 0
SN 0 number of others: 0
SN 0 number of multiallelic sites: 0
SN 0 number of multiallelic SNP sites: 0
因此對于長讀長reads進(jìn)行variant calling操作，還需要更復(fù)雜的工具來做，比如谷歌的DeepVariant 、Ruibang Luo開發(fā)的Clairvoyante
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
