线粒体拥有自身的遗传物质和遗传体系，但其基因组大小有限，是一种半自主细胞器。除了为细胞供能外，线粒体还参与诸如细胞分化、细胞信息传递和细胞凋亡等过程，并拥有调控细胞生长和细胞周期的能力。
电镜下的线粒体.png
线粒体基因.png
来源：
Homo sapiens (ID 51) - Genome - NCBI (nih.gov)
线粒体基因.png
来源：
Genome List - Genome - NCBI (nih.gov)
1.下载线粒体参考序列
wget -c http://ftp.ensembl.org/pub/release-107/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.chromosome.MT.fa.gz
2.构建索引文件
bowtie-build --threads 1 Mus_musculus.GRCm39.dna.chromosome.MT.fa MT
hisat2-build -p 1 Mus_musculus.GRCm39.dna.chromosome.MT.fa MT
参考Hisat2, Bowtie, Bowtie2和BWA构建基因组索引 - 简书 (jianshu.com)
3 过滤
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
4 比对
Bowtie
-n比对模式
ls *1.clean.fq.gz |while read id
do
bowtie -t -p 3 -q -n 3 -l 28 -X 500 -a --best --strata /media/lzx/0000678400004823/Indexs/bowtie/Mm10/mt/MT \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%_*}.bowtie.log\
|samtools sort -@ 2 -o ${id%%.*}_btp.bam
done
-t ，打印每个阶段花费的时间
-p，线程数（或--threads）
-n，-n比对模式，种子区错配的碱基数
-l，左侧至高质量种子区的碱基数，设置越大运行越快
-X，指定fragment的最大长度（最小长度为-I），默认为250
-q，指定输入为fastq文件（-f，指定输入为fasta文件）
-a --best --strata，有多个比对结果时，只报告所有比对结果中匹配度最高的。
-1，指定双端测序中第一个文件
-2，指定双端测序中第二个文件（单端数据用--12）
-S，指定输出sam文件的文件名
报错：
Error while flushing and closing output
terminate called after throwing an instance of 'int'
不知道原因
Hisat2
ls *1.clean.fq.gz |while read id
do
hisat2 -t -q -p 4 \
-x /media/lzx/0000678400004823/Indexs/Hisat2/Mus/MT/MT \
-1 $id -2 ${id%_*}_2.clean.fq.gz 2>${id%_*}.hisat2.log \
|samtools sort -@ 2 -o ${id%_*}_ht2p.bam
done
比对成功。
统计比对结果
samtools flagstat FA51_ht2p.bam
#73199316 + 0 in total (QC-passed reads + QC-failed reads)
#0 + 0 secondary
#0 + 0 supplementary
#0 + 0 duplicates
#66147 + 0 mapped (0.09% : N/A)#共比对上66147条
#73199316 + 0 paired in sequencing
#36599658 + 0 read1
#36599658 + 0 read2
#64870 + 0 properly paired (0.09% : N/A)#成对的有64870条，即32435对。
#65268 + 0 with itself and mate mapped
#879 + 0 singletons (0.00% : N/A)
#0 + 0 with mate mapped to a different chr
#0 + 0 with mate mapped to a different chr (mapQ>=5)
查看比对到染色体的reads
samtools index FA51_ht2p.bam
samtools idxstats FA51_ht2p.bam
#MT 16299 66147 879
#* 0 0 73132290
线粒体长16299 bp，比对上66147条。
samtools fastq -F 4 -1 read1.fq -2 read2.fq FA51_ht2p.bam
samtools view -f 4 -b FA51_ht2p.bam >unmapped.bam
samtools fastq -F 8 -1 read1-1.fq -2 read2-2.fq unmapped.bam
cat read1.fq read1-1.fq >merge1.fq
cat read1.fq read2-2.fq >merge2.fq
seqkit sort -n merge1.fq -o merge1.sort.fq
seqkit sort -n merge2.fq -o merge2.sort.fq
conda install -c bioconda spades
spades.py -1 merge1.sort.fq -2 merge2.sort.fq -o ./spades -t 4 --cov-cutoff auto --careful
