写给女朋友的教程
像人，小鼠这种模式动物，几乎所有比对软件的主页上都可以找到已经构建好的基因组索引文件。要使用的话，直接下载构建好的索引文件就可以了。分析步骤如下：
https://www.jianshu.com/p/eb89ab4af035
https://www.jianshu.com/p/f5e3fab7e9bd
https://www.jianshu.com/p/f9df1bf50b80
但是像一些非模式动物，就需要自己下载和构建索引文件。
NCBI和EMBL和DDBJ都可以下载基因组序列，这次选择去EMBL。
http://asia.ensembl.org/index.html
Fig1.PNG
All genomes，Select a species选择goat，就会跳转到基因组和基因组注释文件下载页面。
Fig2.PNG
左方是基因组，右方是基因组注释文件
点击Download DNA sequence (FASTA)，会进入ftp://ftp.ensembl.org/pub/release-96/fasta/capra_hircus/dna/
ftp://ftp.ensembl.org/pub/release-96/fasta/capra_hircus/dna/，会分别显示出所有染色体的序列文件，点击Capra_hircus.ARS1.dna.toplevel.fa.gz
下载合并好的基因组序列。
但是，校园网下载很慢。
只能用手机流量下，1MB/s
wget ftp://ftp.ensembl.org/pub/release-96/fasta/capra_hircus/dna/Capra_hircus.ARS1.dna.toplevel.fa.gz
NCBI
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/004/361/675/GCA_004361675.1_CVASU_BBG_1.0/GCA_004361675.1_CVASU_BBG_1.0_genomic.fna.gz
点击右方Download GTF，下载基因组注释文件。
基因组索引文件就像一本书的目录，里面集成了比对软件的算法，所以虽然基因组序列是一样的，但不同的比对软件需要不同的基因组索引文件。
基因组注释文件就是把ATCG这些序列按照位置分好区域，告诉你这块区域是什么基因。
构建基因组索引文件：（https://www.jianshu.com/p/89b35626befa）
hisat2-build -p 3 Capra_hircus.ARS1.dna.toplevel.fa Capra_hircus
此处注意-p而非–p，否则回会出现：
Settings:
Output files: "3.*.ht2"
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
Local offset rate: 3 (one in 8)
Local fTable chars: 6
Local sequence length: 57344
Local sequence overlap between two consecutive indexes: 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:
–p
Error: could not open –p
Total time for call to driver() for forward index: 00:00:00
Error: Encountered internal HISAT2 exception (#1)
Command: hisat2-build --wrapper basic-0 –p 3 Capra_hircus.ARS1.dna.toplevel.fa Capra_hircus.ARS1.dna.toplevel
构建成功会提示：
hisat2-build -p 2 Capra_hircus.ARS1.dna.toplevel.fa Capra_hircus.ARS1.dna.toplevel
#Total time for call to driver() for forward index: 00:59:39
生成8个索引文件。
下载高通量测序数据
解压
ls *.sra|while read id;
do
fastq-dump --gzip --split-3 $id;
done
质控
ls *.gz|while read id;
do
fastqc -t 2 $id; #t 线程
done
过滤
mkdir -p fastp
ls *1.fq.gz|while read id;
do
fastp -5 20 -i ${id%_*}_1.fq.gz -I ${id%_*}_2.fq.gz \
-o ${id%_*}_1.clean.fq.gz -O ${id%_*}_2.clean.fq.gz \
-j ./fastp/${id%_*}.json -h ./fastp/${id%_*}.html;
done
比对
ls *1.clean.fq.gz|while read id;
do
hisat2 -t -p 3 \
-x /media/luozhixin/0000678400004823/Indexs/Hisat2/goat/Capra_hircus.ARS1.dna.toplevel \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.hisat2.log \
|samtools sort -@ 2 -o ${id%_*}_ht2p.bam
done
比对结果
Time loading forward index: 00:00:37
Time loading reference: 00:00:04
Multiseed full-index search: 00:16:45
23329781 reads; of these:
23329781 (100.00%) were paired; of these:
772486 (3.31%) aligned concordantly 0 times
21558007 (92.41%) aligned concordantly exactly 1 time
999288 (4.28%) aligned concordantly >1 times
----
772486 pairs aligned concordantly 0 times; of these:
33233 (4.30%) aligned discordantly 1 time
----
739253 pairs aligned 0 times concordantly or discordantly; of these:
1478506 mates make up the pairs; of these:
934201 (63.19%) aligned 0 times
498073 (33.69%) aligned exactly 1 time
46232 (3.13%) aligned >1 times
98.00% overall alignment rate
Time searching: 00:16:49
Overall time: 00:17:26
索引
此步骤是为导入IGV查看所用
只有按染色体位置排序才可以
ls *.bam|while read id;
do
samtools index $id
done
